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The  NSSEFF  research  program  has  been  directed  to  create  new  paradigms  for  the  strong  interaction  of 
light  and  matter.  In  general  terms,  the  research  has  laid  foundations  for  building  complex  quantum  systems 
from  ’simple'  components  of  single  atoms  and  photons.  Theoretical  investigations  have  predicted  surprising 
quantum  phenomena  that  have  not  heretofore  existed  in  Nature.  More  specifically,  the  research 

investigated  new  quantum  phases  of  matter  with  photon-mediated  atom-atom  interactions  in  one  and  twodimensional 
nano-photonic  lattices  (e.g.,  self-organization  of  tree-space  atoms  into  crystals  bound  by  light). 

In  a  complimentary  fashion,  exotic  quantum  phases  for  photons  become  possible  by  way  of  strong  photonphoton 
interactions  created  by  an  underlying  lattice  of  atoms  (e.g.,  two  propagating  photons  bonded  to 
become  a  'molecule'). 

In  terms  of  methodology,  the  NSSEFF  research  required  an  interdisciplinary  'toolkit'  from  atomic  physics, 
quantum  optics,  and  nano-photonics  for  the  control,  manipulation,  and  interaction  of  atoms  and  photons 
with  a  complexity  and  scalability  well  beyond  the  prior  state-of-the-art.  The  NSSEFF  research  achieved 
important  goals  toward  such  integration  across  disciplines  and  utilized  these  tools  to  investigate  novel 
phenomena  in  the  physics  of  light-matter  interactions. 

More  specifically,  with  NSSEFF  support,  a  technical  infrastructure  was  created  for  1)  the  advancement  of 
the  frontier  for  the  design,  fabrication,  and  characterization  of  new  generations  of  1 D  and  2D  photonic 
crystals  in  low-loss  dielectrics,  and  2)  the  integration  of  these  devices  into  the  realm  of  ultra-cold  atomic 
physics  to  produce  nano-scopic  lattices  of  trapped  atoms  that  are  strongly  coupled  to  single  photons  within 
the  photonic  crystals. 

For  operation  in  the  dispersive  regime  of  a  one-dimensional  photonic  crystal  waveguide  (PCW), 
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Atom-light  interactions  in  quasi- ID  nanostructures: 
a  Green’s  function  perspective 
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17  June  2016 

Abstract.  Based  on  a  formalism  that  describes  atom-light  interactions  in  terms 
of  the  classical  electromagnetic  Green’s  function,  we  study  the  optical  response  of 
atoms  and  other  quantum  emitters  coupled  to  one-dimensional  photonic  structures, 
such  as  cavities,  waveguides,  and  photonic  crystals.  We  demonstrate  a  clear  mapping 
between  the  transmission  spectra  and  the  local  Green’s  function  that  allows  to  identify 
signatures  of  dispersive  and  dissipative  interactions  between  atoms,  gaining  insight 
into  recent  experiments. 


PACS  numbers:  42.50.Ct,  42.50.Nn 


Keywords :  Quantum  optics,  nanophotonics,  waveguide  QED. 


1.  Introduction 


As  already  noticed  by  Purcell  in  the  first  half  of  the  past  century,  the  decay  rate  of  an 
atom  can  be  either  diminished  or  enhanced  by  tailoring  its  dielectric  environment  §i 
Likewise,  by  placing  more  than  one  atom  in  the  vicinity  of  photonic  nanostructures,  one 
can  curtail  or  accelerate  their  collective  decay.  In  addition  to  modifying  the  radiative 
decay,  nanophotonic  structures  can  be  employed  to  spatially  and  spectrally  engineer 
atom-light  interactions,  thus  obtaining  fundamentally  different  atom  dynamics  to  those 
observed  in  free-space  |4|. 

In  the  past  decade,  atoms  and  other  quantum  emitters  have  been  interfaced  with  the 
electromagnetic  fields  of  a  plethora  of  quasi- ID  nanostructured  reservoirs,  ranging  from 
high-quality  optical  (h-10|  and  microwave  01  12 1  cavities  to  dielectric  1 13-18] ,  metallic 


f  These  authors  contributed  equally  to  this  research. 
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19  -  22  and  superconducting  1 23 , 24  waveguides.  Photonic  crystal  waveguides,  periodic 


dielectric  structures  that  display  a  bandgap  where  light  propagation  is  forbidden  1 25 , 26 


have  been  proposed  as  promising  candidates  to  study  long-  and  tunable-range  coherent 
interactions  between  quantum  emitters  27-30  .  Due  to  the  different  character  of  the 


guided  modes  at  various  frequencies  within  the  band  structure  of  the  photonic  crystal, 
the  interaction  of  the  quantum  emitters  with  the  nanostructure  can  be  remarkably 
distinct  depending  on  the  emitter  resonance  frequency.  Far  away  from  the  bandgap, 
where  light  propagates,  the  guided  modes  resemble  those  of  a  conventional  waveguide. 
Close  to  the  bandgap,  but  still  in  the  propagating  region,  the  fields  are  similar  to  those 
of  a  quasi- ID  cavity,  whereas  inside  the  bandgap  the  fields  become  evanescent,  decaying 
exponentially. 

All  these  regimes  have  been  recently  explored  in  the  lab,  where  atoms  (31—33  and 
quantum  dots  |4, 34, 35||  have  been  interfaced  with  photonic  crystal  waveguides.  Most 
of  these  experiments  have  been  performed  in  conditions  where  the  resonance  frequency 
of  the  emitter  lies  outside  the  bandgap.  However,  very  recently,  the  first  experiments 
of  atoms  1 36 1  and  superconducting  qubits  [37  interacting  with  evanescent  modes  in  the 
bandgap  of  photonic  crystal  waveguides  have  been  reported. 

Within  this  context,  it  has  become  a  necessity  to  understand  the  rich  spectral 
signatures  of  atom-like  emitters  interacting  through  the  guided  modes  of  quasi  one¬ 
dimensional  nanophotonic  structures  within  a  unified  framework  that  extends  beyond 
those  of  cavity  38  or  waveguide  QED  |39|.  In  this  work,  we  employ  a  formalism  based 
on  the  classical  electromagnetic  Green’s  function  40-44  to  characterize  the  response 
of  atoms  that  interact  by  emitting  and  absorbing  photons  through  the  guided  mode  of 
the  nanostructure.  Since  the  fields  in  the  vicinity  of  the  structure  might  have  complex 
spatial  and  polarization  patterns,  the  full  Green’s  function  is  only  known  analytically  for 


a  handful  of  systems  (such  as  planar  multilayer  stacks  (45 1,  infinite  nanofibers  |46 , 47 


and  a  few  more  |42|)  and  beyond  that  one  has  to  resort  to  numerical  solvers  of  Maxwell’s 
equations.  However,  in  quasi-ID  nanostructures,  one  can  isolate  the  most  relevant 
guided  mode  and  build  a  simple  prescription  for  the  ID  Green’s  function  that  accounts 
for  the  behavior  of  this  mode,  greatly  simplifying  the  problem. 

In  the  first  part  of  the  article,  we  summarize  the  procedure  to  obtain  an  effective 
atom-atom  Hamiltonian,  in  which  the  guided-mode  fields  are  effectively  eliminated  and 
the  atom  interactions  are  written  in  terms  of  Green’s  functions  40-  42 1 .  We  then  apply 


this  formalism  to  a  collection  of  atoms  in  different  quasi  one-dimensional  dielectric 
environments,  and  analyze  the  atomic  transmission  and  reflection  spectra  in  terms 
of  the  eigenvalues  of  the  matrix  consisting  of  the  Green’s  functions  between  every 
pair  of  atoms.  We  show  that,  in  the  linear  (low-saturation)  regime,  asymmetry  in 
the  transmission  spectra  and  frequency  shifts  are  signatures  of  coherent  atom-light 
interactions,  whereas  symmetric  lineshapes  reveal  dissipation.  Finally,  based  on  the 
rapid  technical  advances  in  fabrication  of  both  photonic  and  microwave  structures,  we 
project  observable  signatures  that  can  be  made  in  the  next  generation  of  experiments 
of  atoms  and  superconducting  qubits  interacting  in  the  bandgap  of  photonic  crystal 
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2.  Atom-light  interactions  in  terms  of  Green’s  functions 

Much  effort  has  gone  into  developing  a  quantum  formalism  to  describe  atoms  coupled  to 
radiation.  A  conventional  technique  is  to  express  the  field  in  terms  of  a  set  of  eigenmodes 
of  the  system,  with  corresponding  creation  and  annihilation  operators  af  and  a  [39].  This 
canonical  quantization  technique  is  well  suited  for  approximately  closed  systems  such 
as  high-Q  cavities  and  homogeneous  structures  such  as  waveguides,  both  of  which  have 
simple  eigenmode  decompositions.  However,  the  application  of  this  quantization  scheme 
to  more  involved  nanostructures  is  not  straight  forward.  Further,  the  formalism  is  not 
suited  for  dispersive  and  absorbing  media  as  the  commutation  relations  for  the  field 
operators  are  not  conserved  [48|. 

Instead,  here  we  describe  atom- light  interactions  using  a  quantization  scheme  based 
on  the  classical  electromagnetic  Green’s  function,  valid  for  any  medium  characterized  by 
a  linear  and  isotropic  dielectric  function  e(r,  uj),  closely  following  the  work  of  Welsch  and 
colleagues  140-42,44  .  In  the  following,  we  employ  this  formalism  to  derive  an  atom- 


atom  Hamiltonian  in  which  the  field  is  effectively  eliminated,  yielding  an  expression 
that  only  depends  on  atomic  operators.  Moreover,  once  the  dynamics  of  the  atoms  is 
solved,  the  electric  field  at  every  point  along  the  quasi  one-dimensional  structure  can  be 
recovered  through  an  expression  that  relates  the  field  to  the  atomic  operators. 

Classically,  the  field  E(r,  uj)  at  a  point  r  due  to  a  source  current  j(r',o;)  at  r'  is 
obtained  by  means  of  the  propagator  of  the  electromagnetic  field,  the  dyadic  Green’s 
function  (or  Green’s  tensor),  as  E(r,  uj)  =  i fj,0oj  f  dr'  G(r,  r',  uj)  ■  j(r',o;).  In  particular, 
for  a  dipole  source  p  located  at  r0,  the  current  is  j(r,o;)  =  — kupfi(r  —  r0),  and  the  field 
reads  E(r,a;)  =  hqoj2  G(r,  r0,  uj)  ■  p.  The  tensorial  structure  of  the  Green’s  function 
accounts  for  the  vectorial  nature  of  the  electromagnetic  field,  as  a  dipole  directed  along 
the  ^-direction  can  create  a  field  polarized  not  only  along  x,  but  also  along  y  and  z  [§J 

The  Green’s  function  G(r,  r',  oj)  is  the  fundamental  solution  of  the  electromagnetic 
wave  equation,  and  obeys  [49]: 


u 


V  x  V  x  G(r,  r',  oj) - -e(r,  oj)  G(r,  r',  oj)  =  S(r  —  r')l, 


(1) 


where  e(r ,oj)  is  the  medium  relative  permittivity.  For  a  scalar  permittivity,  Lorentz 
reciprocity  holds  and,  then,  GT(r,  r',oj)  =  G(r',r ,oj),  where  T  stands  for  transpose 
(and  operates  on  the  polarization  indexes).  In  analogy  to  its  classical  counterpart, 
the  electric  field  operator  can  be  written  in  terms  of  bosonic  annihilation  (creation) 
operators  f  (f^)  as 


40 


E(r,cn)  =  i/i0  uj2J^^-  J  dr1  y/lm^r',  cu)}  G(r,  r'oj)  •  f^w)  +  h.c. 


(2) 


—  E~*~(r, oj)  +  E  (r,cn), 


Throughout  this  manuscript., the, fimen’.s. tensor  will  be  also,  denoted  as  Green’s  function. 
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where  E+(_)(r,a;)  is  the  positive  (negative)  frequency  component  of  the  field  operator, 
and  h.c.  stands  for  Hermitian  conjugate.  Within  this  quantization  framework, 
f(r,  uj)  is  associated  with  the  degrees  of  freedom  of  local  material  polarization  noise, 
which  accompanies  the  material  dissipation  Im{e(r,  eu)}  as  required  by  the  fluctuation- 
dissipation  theorem  44  .  This  expression  guarantees  the  fulfillment  of  the  canonical 
field  commutation  relations,  even  in  the  presence  of  material  loss.  The  appearance  of 
the  Green’s  function  reveals  that  the  quantumness  of  the  system  is  encoded  in  either  the 
correlations  of  the  noise  operators  f  or  in  any  other  quantum  sources  (such  as  atoms), 
but  the  field  propagation  obeys  the  wave  equation  and  as  such  the  spatial  profile  of  the 
photons  is  determined  by  the  classical  propagator. 

We  now  want  to  investigate  the  evolution  of  N  identical  two-level  atoms  of 
resonance  frequency  oj_\  that  interact  through  a  guided  mode  probe  field  of  frequency  u>p. 
Within  the  Born-Markov  approximation,  we  trace  out  the  photonic  degrees  of  freedom, 
obtaining  an  effective  atom- atom  Hamiltonian  41,50,51  .  This  approximation  is  valid 


when  the  atomic  correlations  decay  much  slower  than  the  photon  bath  correlations, 
or,  in  other  words,  when  the  Green’s  function  is  characterized  by  a  broad  spectrum, 
which  can  be  considered  to  be  flat  over  the  atomic  linewidth.  Then,  the  atomic  density 
matrix  evolves  according  to  px  =  — (i /h)  [5£,pa]  +  2?[pa]  |38|.  Within  the  rotating 
wave  approximation,  and  in  the  frame  rotating  with  the  probe  field  frequency,  the 
Hamiltonian  and  Lindblad  operators  read 


*  =  -fiA*  £  4  -  ft  £  -£(d-  E“(r,)  4  +  d-  •  E+(r()  a%)  , 

2=1  2,J  =  1 


N 


N 

E 

2=1 


l  J  /  . 

&[Pa\  =  I]  -it  {2algepAaJeg  -  alegaJgepA  -  pAoleQV 


eg  ge 


i,j= 1 


). 


(3a) 

(3b) 


where  Ep  is  the  guided  mode  probe  field,  and  A  a  =  u>p  —  u>a  is  the  detuning  between 
the  guided  mode  probe  field  and  the  atom.  The  dipole  moment  operator  is  expressed 
in  terms  of  the  dipole  matrix  elements  as  p;  =  d *  <rJeg  +  d  aJge ,  where  a:lg  =  \e)  (g\  is 
the  atomic  coherence  operator  between  the  ground  and  excited  states  of  atom  j ,  and 
d  =  (g|d|e)  is  the  dipole  matrix  element  associated  with  that  transition.  The  spin- 
exchange  and  decay  rates  are 


Jv  =  (po^p/fi)  d*  •  Re  G(ri;  r j,ujp)  ■  d,  (4a) 

ry  =  (2p0  uip/h)  d*  •  ImG^i.r^Wp)  •  d.  (4b) 

Note  that  the  dispersive  and  dissipative  atom-atom  couplings  are  given  in  terms  of 
the  total  Green’s  function  of  the  medium.  For  a  quasi-ID  nanostructure,  the  Green’s 
function  can  be  expressed  as  G(rj,  iq, u>p)  =  GiD(rj,  iq, u>p)  +  G'(rj,  iq,  cnp),  where  the 
first  term  corresponds  to  the  guided  mode  that  propagates  along  the  structure  mediating 
atom-atom  interactions,  and  the  second  term  accounts  for  all  other  field  modes  (e.g., 
emission  into  free  space). 
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Due  to  the  fast  spatial  decay  of  the  non-guided  Green’s  function,  the  interaction 
mediated  by  G'(rj,  i*j,  o;p)  is  not  collective  in  the  low-density  limit  (i.e.  when  the 
atoms  are  far  away  from  each  other).  We  can  then  write  J'J  =  +  J'St]  and 

pi  _  p^D  _|_  T'Sij,  where  d,j  is  the  Kronecker  delta.  In  particular,  in  free-space,  F' 
is  simply  r0  =  (2/i0Wp/h)d*  •  Im  G0(rj,  r*,  cnp)  •  d  =  (Up|d|2/37rhe0c3,  where  G0  is  the 
vacuum’s  Green’s  function  [i.e.  the  solution  to  Eq.  ([!])  when  e(r,  u)  =  1].  Depending  on 
the  geometry  and  dielectric  response  of  the  nanostructure,  and  on  the  atom  position,  T' 
can  be  larger  or  smaller  than  r0.  J’  accounts  for  frequency  shifts  due  to  other  guided 
and  non-guided  modes,  and  is  in  general  spatially  dependent.  In  fact,  the  value  of  J'  is 
dependent  upon  particular  details  of  atom  trapping  and  geometry  of  the  nanostructure. 
We  will  for  simplicity  consider  J'  identical  for  every  atom  and  assume  that  this  constant 
value  has  been  incorporated  into  the  definition  of  uja- 

Once  the  dynamics  of  the  atomic  coherences  are  solved  for,  one  can  reconstruct  the 
field  at  any  point  in  space.  Generalizing  Eq.  (6.16)  of  Ref.  |4l)  for  more  than  a  single 
atom,  the  evolution  of  the  bosonic  field  operator  is  given  by 


N 


f(r,w)  =  -iwf(r,w)  +  ^  ^ -J—Im{e(r,  uj)}  ^  G*(r,  Tj,u)  •  d  aJge, 


(5) 


where  the  atoms  act  as  sources  for  the  bosonic  fields.  We  can  formally  integrate  this 
expression  and  plug  it  into  the  equation  for  the  field  [Eq.  ([2])].  After  some  algebra, 
and  performing  Markov’s  approximation,  we  arrive  at  the  final  expression  for  the  field 
operator,  which  is  simply 


N 


E+(r)  =  E+(r)+«,^5:G(r,rJ, 

3  = 1 


(6) 


This  expression  can  be  understood  as  a  generalized  input-output  equation,  where  the 
total  guided  mode  field  is  the  sum  of  the  probe,  i.e.  free,  field  Ep  (r)  and  the  field  re¬ 
scattered  by  the  atoms.  The  quantum  nature  of  these  equations  has  been  treated  before 
when  deriving  a  generalized  input-output  formalism  for  unstructured  waveguides  [52 , 53 1 . 


3.  Transmission  and  reflection  in  quasi-ID  systems 


3.1.  Atomic  coherences  in  the  low  saturation  regime 

We  now  explore  the  behavior  of  the  atoms  under  a  coherent,  continuous- wave  probe  field. 
In  the  single-excitation  manifold  and  low  saturation  (linear)  regime  ((cree)  =  0),  the 
atoms  behave  as  classical  dipoles.  Then,  the  Heisenberg  equations  for  the  expectation 
value  of  the  atomic  coherences  ({&eg)  =  creg)  are  linear  on  the  atomic  operators,  and 
read 


i  9ij 


rr 


gel 


where  D*  =  d*  •  E+  (r,;) / h  is  the  guided  mode  Rabi  frequency  (with  Ep  =  (Ep)),  and 
gij  =  AT  ' d 


(7) 
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depends  only  on  the  Green’s  function  of  the  guided  mode.  For  long  times,  the  coherences 
will  damp  out  to  a  steady  state  ( &lge  =  0).  The  solution  for  the  atomic  coherences  is 
then 


age 


i  i 


SI  with  M  =  (Aa  +  ir'/2)  1  +  0. 


(8) 


In  the  above  equation,  age  =  (cde,  ...  ,cr^e)  and  O  =  (Qi,  ...  ,Qn)  are  vectors  of 
N  components,  and  JH  is  a  N  x  N  matrix  that  includes  the  dipole-projected  matrix 
0  of  elements  <7^.  Significantly,  the  matrix  is  not  Hermitian,  as  there  is  radiation 
loss.  However,  due  to  reciprocity,  the  Green’s  function  matrix  is  complex  symmetric 
[GT(r,  y',uj)  =  G(r',  r.  u;)],  and  0  inherits  this  property  if  the  dipole  matrix  elements 
are  real,  which  will  be  a  condition  enforced  from  now  on.  Complex  symmetric  matrices 
can  be  diagonalized,  0v^  =  with  £  =  1 . . .  N,  where  A^  and  are  the  eigenvalues 
and  eigenvectors  of  0,  respectively.  Since  the  first  term  of  M  is  proportional  to  the 
identity,  Jli  and  0  share  the  same  set  of  eigenvectors. 

The  eigenmodes  represent  the  spatial  profile  of  the  collective  atomic  excitation, 
i.e.,  the  dipole  amplitude  and  phase  at  each  atom.  However,  as  the  matrix  0  is  non- 
Hermitian,  the  eigenmodes  are  not  orthonormal  in  the  regular  sense,  but  instead  follow 
different  orthogonality  and  completeness  prescriptions,  namely  vj  ■  v^/  =  8^^  and 


£  V=1 


rT  - 


3L,  where  T  indicates  transpose  instead  of  the  customary  conjugate 


transpose  54  .  After  inserting  the  completeness  relation  into  Eq.  (|8h ,  we  find  that  the 


expected  value  of  the  atomic  coherences  in  the  steady  state  in  terms  of  the  eigenvalues 
and  eigenvectors  of  the  quasi-ID  Green’s  function  is 


&ge  ~ 


-  E  7T- 

(Gmode  V  A 


(vf 


•  n) 


-h 


^,id)  +  i(p  +  n,iD)/2 


(9) 


where  =  He  A^  and  T^id  =  2  Im  A^  are  the  frequency  shifts  and  decay  rates 

corresponding  to  mode  £,  and  the  sum  is  performed  over  mode  number  from  1  to  N.  The 
scalar  product  in  the  numerator  vj  •  Q  =  J2f=  1  vz,j  %  describes  the  coupling  between 
the  probe  field  and  a  particular  collective  atomic  mode. 

Therefore,  the  dynamics  of  the  atoms  can  be  understood  in  terms  of  the  eigenmodes 
of  0,  where  the  real  and  imaginary  parts  of  the  eigenvalues  correspond  to  cooperative 
frequency  shifts  and  decay  rates  of  the  collective  atomic  modes  {£}.  As  the  modes  are 
non-normal,  the  observables  cannot  be  expressed  as  the  sum  over  all  different  mode 
contributions  but,  instead,  any  measurable  quantity  will  show  signatures  of  interference 
between  different  modes.  Although  it  could  be  considered  a  mathematical  detail,  the  fact 
that  the  modes  of  a  system  are  non-normal  has  deep  physical  consequences.  For  instance, 
non-normal  dynamics  is  responsible  of  phenomena  as  different  as  the  Petermann  excess- 


noise  factor  observed  in  lasers  {55-  57  or  the  transient  growth  of  the  shaking  of  a  building 
after  an  earthquake  58  . 
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3.2.  Transmission  and  reflection  coefficients 


Having  previously  calculated  the  linear  response  of  an  ensemble  of  atoms  to  an  input 
field,  we  now  relate  the  response  to  observable  outputs,  i.e.  the  reflected  and  transmitted 
fields.  One  can  calculate  the  total  field  from  Eq.  by  substituting  in  the  solution  of 
Eq.  ([£])  for  the  atomic  coherences  age.  For  the  sake  of  simplicity,  we  now  assume  that 
the  atomic  chain  and  the  main  axis  of  the  nanostructure  are  oriented  along  x,  and  the 
atoms  have  all  the  same  radial  position  and  thus  the  same  ’transversal’  coupling  into 
the  nanostructure.  The  field  is  considered  to  be  polarized  along  y,  and  reads 


Vf  -K 


N  (gT(x)  •  A ,T 

E+(x )  =  Efix)  =  Eflix)  —  — — - r — — - — — , 

v  P  fri(  Aa  +  4id)  +  i(r'  +  rSjlD)/2 


(10) 


where  the  j-component  of  vector  g(x)  is  gflx)  =  g(x,Xj)  =  (/xocOpd2 /h)Gu),yy(x,Xj,cup), 
and  gT  (x)  •  =  Y^Lj  9j{x)v£,j  represents  how  much  the  mode  £  contributes  to  the  field 
emitted  by  the  atoms.  Note  that  now  the  electric  field  vector  E+  in  Eq.  (10)  no  longer 
represents  different  polarization  components,  but  a  single  polarization  at  different  atom 
positions. 

In  order  to  connect  the  above  expression  to  the  transmission  and  reflection 
coefficients,  we  evaluate  the  field  E+(x)  at  the  positions  x  =  :i'rigilt  and  x  =  uqeft  ,  which 
are  considered  to  be  immediately  outside  the  atomic  chain.  The  details  of  the  derivation 
are  provided  in|Appendix  A|  The  normalized  transmission  and  reflection  coefficients  are 


*(Aa)Ao(Aa)  =  1 


1 

g (bright)  Teft ) 


N 


E 

£=1 


(gT(^right)  '  V€)  (vj  ■  g(xief t)) 

(^a  +  ^,id)  +  +  r^m)/2 


KAa)  =  r„(AA)  _  1  y  •  y)  (v[ ■  eM) 

^(^left^left)  g=1  (Aa  +  T^id)  +  i(T' +  T^;id)/2’ 


(11a) 

(lib) 


where  t0(A\)  and  t0(Aa)  are  the  transmission  and  reflection  coefficients  for  the  ID 
photonic  structure  when  no  atoms  are  present.  One  can  further  simplify  the  expression 
for  the  transmission  so  that  the  resulting  equation  only  depends  on  the  eigenvalues,  and 
not  the  eigenfunctions,  of  the  Green’s  function  matrix  0  (as  shown  in  Appendix  B). 
Then, 


Appendix  B 


N 


t(AA)/t0(AA)  =  JJ 
(=1 


_ Aa  +  ir/2 _ 

(Aa  +  Aid)  +  i(r;  +  Aid)/2 


N 


n  ^(Aa)- 

€=1 


(12) 


The  total  transmission  coefficient  can  thus  be  written  as  the  product  of  the  transmission 
coefficients  of  each  of  the  collective  atomic  modes.  Noticeably,  when  looking  at  the 
transmission  spectrum  of  atoms  that  interact  through  the  guided  mode  of  a  quasi-ID 
nanostructure,  there  is  a  redundancy  between  the  eigenfunctions  and  eigenvalues,  and 
one  is  able  to  obtain  an  expression  that  does  not  depend  on  the  former  (i.e.  all  the 

relevant  information  about  the  geometry  is  contained  in  the  collective  frequency  shifts 
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(a) 


c  •  J’’ 


F 


Figure  1.  (a)  Sketch  of  an  atom  interacting  with  the  guided  mode  of  a  structured 
ID  nanostructure.  The  single-atom  decay  rate  is  T id,  and  the  decay  into  non-guided 
modes  is  characterized  by  T'.  (b)  Normalized  transmission  spectra  (|t/to|2)  for  a  single 
atom  for  different  values  of  the  ratio  between  the  real  and  imaginary  parts  of  the  guided 
mode  Green’s  function,  following  Eq.  (13).  The  decay  rate  into  the  guided  modes  is 
taken  to  be  Tid  =  T'  for  all  cases. 


and  decay  rates).  In  particular,  for  a  single  atom  located  at  Xj  with  =  Jid  and 
r%  =  r1D,  the  eigenvalues  are  directly  proportional  to  the  local  Green’s  function,  and 


£(Aa)/£0(Aa)  = 


Aa  +  lT/2 


(Aa  +  Jid)  +  i(F  +  Tid)/2 


(13) 


The  transmittance  T  =  \t\2  can  be  recast  into  a  Fano-like  lineshape  59 

2 


as 


r/To  =  + 

l  +  x2 


r 


i 


r'  +  r1D;  i  +  x2 


(14) 


where  x  =  2(Aa  +  JiD)/(riD  +  r/)  and  q  =  —  2  Jid/ (Fid +  F)  is  the  so-called  asymmetry 
parameter.  For  T'  <C  Tid,  the  second  term  is  negligible  and  T/T0  is  a  pure  Fano 
resonance,  with  q  =  —  Re { G  i  d  ( r?  •  r?  •  u;p ) } / Im { G  i d  ( iq ,  r, ,  cjp ) }  ■  Fano  resonances  arise 
whenever  there  is  interference  between  two  different  transport  channels.  For  instance, 
in  a  cavity  far  from  resonance,  there  is  interference  arising  from  all  the  possible  optical 
paths  that  contribute  to  the  transmission  signal  due  to  reflections  at  the  mirrors,  whereas 
in  an  unstructured  waveguide  there  is  no  such  interference  and  thus  the  lineshape  is 
Lorentzian. 

For  a  single  atom,  there  is  a  clear  mapping  between  the  spectrum  lineshape  and 
the  local  ID  Green’s  function.  For  a  nanostructure  with  a  purely  imaginary  self  Green’s 
function  G {xi,xf)  (such  as  a  wave- guide  or  a  cavity  at  resonance),  the  spectrum  is 
Lorentzian,  and  centered  around  the  atomic  frequency.  However,  if  the  real  part  is 
finite,  one  would  observe  a  frequency  shift  of  the  spectrum,  which  becomes  asymmetric. 
Figure  [l](b)  shows  how  the  normalized  transmission  spectrum  for  a  single  atom  becomes 
more  and  more  asymmetric  for  higher  ratios  Jid /Ido •  Also,  there  is  an  appreciable 

blueshift  of  the  spectral  features.  _  .  x„  . 
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We  would  like  to  remark  that  the  Markov  approximation  has  thus  far  been  employed 
in  our  analysis,  as  every  Green’s  function  is  considered  to  be  a  complex  constant  over 
frequency  ranges  larger  than  the  linewidth  of  the  atoms.  If  that  is  not  the  case,  it  is 
not  possible  to  find  simple  expressions  for  the  Hamiltonian  and  Lindblad  terms  for  the 
atomic  density  matrix.  However,  the  expressions  for  the  transmission  and  reflection 
coefficients  are  valid  even  when  the  spectral  variation  of  the  Green’s  function  occurs 
within  frequency  intervals  comparable  to  and  smaller  than  the  atomic  linewidth.  This 
fact  might  not  be  surprising  as,  in  the  low  saturation  limit,  atoms  behave  as  classical 
dipoles,  and  an  equation  for  the  transmission  coefficient  identical  to  Eq.  (12)  can  be 
found  for  classical  emitters,  without  resorting  to  Markov’s  approximation. 


4.  Application  to  several  one-dimensional  photonic  structures 

In  this  section,  we  analyze  the  transmission  spectra  of  atoms  placed  along  common 
quasi-ID  nanostructures,  such  as  cavities,  waveguides,  and  photonic  crystals. 


4-1.  Standing-wave  cavities 

To  begin  with,  we  want  to  illustrate  the  connection  between  the  Green’s  function 


formalism  and  the  well-known  Jaynes  Cummings  (JC)  model  38  .  For  N  atoms  in  a 


driven  cavity  of  length  L  and  effective  area  A,  the  JC  Hamiltonian,  and  its  corresponding 
Lindblad  operator  read 


N  N 

=  —hAcala  -  HAa  +  (atd*e  +  (tlega)  +  hp  (a  +  a1'), 


1=1 


1=1 


■v  N 


2[/3]  =  —  £  {2aleP°Jep  ~ 

z  a*=i 


ger  eg  ~  eg~  ge 


\eP  ~  P&'eg&ge)  +  y  (^apo)  ~  of  Op  ~  pafa)  , 


*0= 


(15a) 

(15b) 


where  a  is  the  cavity-field  annihilation  operator,  p  is  the  density  matrix  for  the  atoms 
and  the  cavity  field,  p  is  a  frequency  that  represents  the  amplitude  of  the  classical  driving 
field,  Ac  =  (Up  —  ojc  is  the  detuning  between  the  driving  (probe)  and  the  cavity  fields, 
and  nc  is  the  cavity- field  decay.  The  atom  cavity  coupling  is  <pt  =  cos(A:c.i:(),  where 

<p  =  d^ujc/ (heoLA)  is  modulated  by  a  function  that  depends  on  the  atoms’  positions 
and  the  cavity  wave- vector  kc.  The  Heisenberg  equations  of  motion  for  the  field  and 
atomic  operators  are 


N 

«  -  i  H  <hi°lge 
1=1 


~  lTh 


a, 


ge 
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When  T'  <C  nc  and  <  min{Ac, /-cc},  the  cavity  held  can  be  adiabatically  eliminated, 
and  the  held  operator  re-expressed  in  terms  of  the  atomic  ones,  i.e. , 


<  i  —  0  — y  a  — 


(Ac  +  i^) 


(  N 

V  +  Y  W'g 

\  i=l 


9e 


Introducing  this  expression  back  into  the  equation  for  the  atomic  operator,  one  can 
deduce  a  master  equation  for  the  atomic  density  matrix  p\.  The  new  Hamiltonian 
and  Lindblad  operators  read  just  as  those  of  Eqs.  (3a)  and  (3b),  but  for  a  classical 
driving  held,  and  with  spin  exchange  and  decay  rates  into  the  cavity  mode  given 
by  JIdOTd)  =  Jm(rw)cos(kcXi)cos(kcXj),  with  Jm  =  — q,2Ac/(A2  +  k2JA),  and 
Tid  =  <p2nc/{ A2  +  k2/4).  It  can  thus  be  seen  that  the  Markovian  approximation  to 
arrive  at  these  equations  is  equivalent  to  the  absence  of  strong  coupling  effects  within 
the  JC  model. 

The  last  step  for  connecting  this  simple  model  with  our  formalism  is  to  calculate  the 
Green’s  function  of  a  cavity  and  conhrm  that  Jff  and  TjJD  are  precisely  those  obtained 
within  the  JC  framework.  The  Green’s  function  of  a  quasi- ID  cavity  formed  by  partially 
transmitting  mirrors  of  reflection  coefficient  r  (chosen  to  be  real)  is  |60| 


CuD  (*Ei  j  3Cj  i 


jc2  eikp\xi-Xj\  _|_  reikp(L+Xi+Xj)  _|_  reikp[L-(xi+Xj)]  _|_  r2^ikp{2L-\xi-Xj\) 


2vgUJpA 


/^-»2g2iA/pX/ 


where  vg  is  the  group  velocity.  For  high-Q  standing- wave  cavities,  i.e.  with  r  ~  1,  and 
choosing  vg  =  c,  the  Green’s  function  can  be  approximated  as 


GiD(xi,Xj,cjp)  ~ 


2ic 


1 


k(jJpA  )  1  —  r2e2ikPL 


cos (kpXi)  cos (kpXj). 


The  cavity  is  resonant  at  a  frequency  u>c  with  corresponding  wave- vector  kc ,  chosen  to 
be  such  that  kcL  =  2ir m,  with  m  being  an  integer.  Close  to  resonance,  one  can  write 
kp  =  kc  +  fik,  and  assume  that  8kL  <C  1.  Then  1  —  r2e2lkpL  ~  1  —  r2  —  2ir28kL,  and  the 
Green’s  function  is  simply 


Xj ,  £Up)  — 


1 


^ujpLA  J  Ac  -)-  1kc/2 


cos (kcXi)  cos (kcXj), 


where  «c  =  (1  —  r2)c/L  is  the  cavity  linewidth.  Therefore,  the  atoms’  spin-exchange 
and  decay  rates  are  given  by 


•/id  =  (po^pd  /h)  ReGw(xi,Xj,ujp)  = 


Tfo  =  (2/j,0uj2d2 /ti)  Im  Gm(xi,  Xj,u>p)  = 


Ar 


(A2  +  k2/4) 


'yi'r3  (A2  +  k?/4) 

which  is  precisely  what  is  obtained  within  the  Jaynes  Cummings  model. 


=  Jid  cos(kcXi)  cos(/ccXj), 


=  Tid  cos (kcXi)  cos (kcXj), 
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Let’s  now  look  at  the  transmission  spectrum  of  N  atoms  in  a  cavity.  As  we  have 
just  demonstrated,  coefficients  of  the  dipole-projected  Green’s  function  matrix  g  read 


gij  =  g(u p)  cos (kcXi)  cos (kcXj), 


(17) 


where  g(u>p)  =  Jib  +  iEio/2.  Depending  on  the  detuning  between  the  probe  field  and 
the  cavity  resonance,  g(ujp)  can  be  purely  imaginary,  yielding  dissipative  atom-atom 
interactions,  or  can  have  both  real  and  imaginary  parts,  resulting  in  both  dissipative 
and  dispersive  couplings. 

The  matrix  g  is  separable  (has  rank  one)  as  it  can  be  written  as  the  tensor  product 
of  just  one  vector  by  itself.  The  matrix  has  one  eigenstate  describing  a  superposition 
of  atomic  coherences  that  couples  to  the  cavity  (a  "bright  mode"),  with  eigenvalue 
Ab  =  YliLi  gn  =  (Jid  +  iTic/2)  Y,f=i  cos 2{kcXi).  This  atomic  collective  excitation  follows 
spatially  the  mode  profile  of  the  cavity,  i.e.  alge  oc  cos (kcXi).  The  matrix  g  has  also  N—  1 
decoupled  ("dark")  modes  of  eigenvalue  0.  Because  these  dark  modes  have  a  zero  decay 
rate  into  the  cavity  mode,  it  is  also  impossible  to  excite  them  employing  the  cavity  field. 
The  optical  response  is  thus  entirely  controlled  by  the  bright  mode,  and  the  transmission 
is  simply 


t(AA)/t0(AA) 


Aa  +  iT'/2 

(Aa  +  Eg !  Jib)  +  i(r  +  Egi  r«D)/2  • 


(18) 


Remarkably,  this  expression  is  valid  no  matter  the  separation  between  the  atoms 
or  whether  they  form  an  ordered  or  disordered  chain.  The  transmission  spectrum 
corresponds  to  that  of  a  ‘super- atom’,  where  the  decay  rates  and  the  frequency  shifts 
are  enhanced  (N-fold  if  all  the  diagonal  components  of  g  are  equal)  compared  to  those  of 
a  single  atom.  This  result  replicates  the  well-known  expressions  for  conventional  cavity 
QED. 


4-2.  Waveguides 


Another  paradigm  that  has  been  investigated  frequently  is  that  of  "waveguide  QED"  39 


The  simple  model  of  such  a  system  consists  of  a  single  guided  mode  with  translational 
invariance,  and  where  the  dispersion  relation  is  well-approximated  as  linear  around  the 
atomic  resonance  frequency.  In  a  ID  translationally  invariant  system,  a  source  simply 
emits  a  plane  wave  whose  phase  at  the  detection  point  is  proportional  to  the  distance 
of  separation.  Therefore,  the  elements  of  the  Green’s  function  matrix  g  depend  on  the 
distance  between  the  atoms,  and  read 


9ii  = 


(19) 


Remarkably,  the  self  Green’s  function  in  a  waveguide  is  purely  imaginary.  The  coherent 
interactions  between  atom  i  and  atom  j  are  dictated  by  the  Hamiltonian  [given  by 
Eq.  (3a)],  and  are  proportional  to  Rejgq}  =  —  (T /2)  sin  kp\xi  —  Xj |,  whereas  the 
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Figure  2.  (a)  Frequency  shifts  and  (b)  decay  rates  of  the  collective  modes  of  a  regular 
chain  of  5  atoms  placed  along  a  waveguide  normalized  to  the  single-atom  decay  rate 
into  the  guided  mode  r id,  as  a  function  of  the  distance  d  between  the  atoms  in  units 
of  the  probe  wavelength. 


dissipation  is  given  by  the  Lindblad  operator  [given  by  Eq.  (3b)],  which  is  proportional 
to  Im {gij}  =  (Tid/2)  coskp(xi  —  xf)  1 22 , 50  .  It  is  thus  clear  that  by  carefully  tuning 


the  distance  between  the  emitters,  one  can  engineer  fully  dissipative  interactions.  If  the 
atoms  form  a  regular  chain  and  are  spaced  by  a  distance  d  such  that  kpd  =  mr,  where 
n  is  an  integer  number,  the  matrix  g  has  only  one  non-zero  eigenvalue  Ab  =  iArr1£)/2 
associated  with  the  bright  atomic  mode.  This  situation  is  analogous  to  the  case  of 
atoms  interacting  in  an  on- resonance  cavity.  Therefore,  there  will  not  be  any  collective 
frequency  shift,  and  the  lineshape  will  be  a  Lorentzian  of  width  Tb  +  T'.  For  n  even, 
the  phases  of  the  dipole  moments  of  the  atoms  are  all  identical,  whereas  for  odd  n  the 
dipole  moments  of  adjacent  atoms  are  7r  out  of  phase. 

For  a  regular  chain  with  lattice  constant  different  from  kpd  =  rm,  or  for  atoms 

placed  randomly  along  the  waveguide,  the  coefficients  of  matrix  g  have  both  a  real 

and  imaginary  part,  and,  to  the  best  of  our  knowledge,  there  is  no  analytic  expression 

for  the  eigenvalues  of  g.  Figure  [2]  shows  the  frequency  shifts  and  decay  rates  of  the 

collective  modes  of  a  N  =  5  atom  chain  as  a  function  of  the  separation  between  the 

atoms.  For  separations  where  kpd  =  rm,  the  real  part  of  the  Green’s  function  is  zero 

and  the  imaginary  part  of  all  modes  but  one  goes  to  zero,  whereas  for  other  spacings  one 

generically  gets  a  zoo  of  coherent  and  dissipative  couplings  of  comparable  strength.  This 

occurs  because  the  real  and  imaginary  parts  of  gl3  are  generically  of  similar  magnitude. 

Figure  [3]  shows  the  transmission  and  reflection  spectra  for  N  =  20  atoms  separated  by 
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Figure  3.  (a)  Normalized  transmission  spectra  for  20  atoms  interacting  through 

the  guided  modes  of  an  unstructured  waveguide.  The  blue  line  represents  a  regular 
separation  between  the  atoms  of  d  =  Ap/2.  The  orange  curves  show  10  different  spectra 
obtained  by  randomly  placing  the  atoms  along  the  nanostructure.  The  black  curve 
represents  the  "non-interacting"  case  of  Eq.  (J20|.  (b)  Normalized  reflection  spectra  for 
the  same  situations  as  in  (a).  We  have  chosen  Tid  =  T'. 


kpd  —  7 r  (blue  curve),  and  for  several  random  realizations  where  each  atomic  position 
is  chosen  randomly  from  a  distribution  kpXi  E  [0,  2tt]  (orange  curves).  The  black  line 
represents  the  non-interacting  case,  which  is  obtained  by  setting  the  non-diagonal  terms 
of  0  to  zero,  yielding  a  transmission  spectrum 

\  N 


£(Aa)Ao(Aa)  = 


Aa  +  ir'/2 


A, 


(20) 


i(r  +  r1D)/2y  ’ 

where  the  transmission  coefficient  is  a  product  of  the  transmission  coefficient  of  each 
singie  atom,  and  the  frequency  shifts  and  decay  rates  are  not  coiiective  quantities  but, 
instead,  singie- atom  parameters. 

Figure  [3](  a)  also  shows  that,  for  random  filling,  although  the  atoms  interact  with 
each  other  (gyy,t  ^  0),  the  transmission  spectra  follow  closely  that  of  a  non-interacting 
system,  for  which  all  the  off-diagonal  elements  are  zero  (gyy*  =  0),  and  the  eigenvalues  of 
matrix  g  are  proportional  to  the  self  Green’s  functions  [G(.i:,:,  x*)]  at  the  atoms’  positions. 
In  this  case,  the  behavior  of  the  emitters  cannot  be  understood  in  terms  of  the  ‘super¬ 
atom’  picture,  as  the  transmission  spectrum  of  the  system  is  significantly  different  from 
a  Lorentzian.  In  particular,  for  the  non-interacting  scenario,  one  can  recast  Eq.  (20)  into 
an  exponential,  and  the  transmittance  recovers  the  well-known  form  of  a  Beer-Lambert 
law,  reading 


r(AA)/r0(AA)  =  exp 


„,nAi  +  <r'  +  riD)74l 

OD 

Ai  +  F2/4 

—  LXjJ 

1  +  (2Aa/F)2_ 

(21) 


where  OD  =  2AVid/T'  is  the  optical  depth  and  the  last  equality  holds  for  Tid  "C  T'. 
This  is  exactly  the  same  behavior  that  an  atomic  ensemble  in  free  space  would  exhibit. 
This  occurs  only  for  non- negligible  T',  which  suppresses  multiple  reflections.  Otherwise 
one  would  see  huge  fluctuations,  .associated  with  Anderson  localization  in  the  spectra. 
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Figure  4.  Collective  frequency  shifts  of  the  modes  of  a  regular  chain  of  N=10  atoms 
in  the  bandgap  of  an  infinite  photonic  crystal  as  a  function  of  nxd,  where  k~  1  is  the 
spatial  range  of  the  interaction  and  d  is  the  distance  between  atoms.  The  atoms  are 
placed  at  even  antinodes  of  the  Bloch  modes. 


The  reflectance  spectrum,  on  the  other  hand,  is  more  complex  and  carries  more 
information  than  the  transmittance,  as  shown  in  Fig.  |3|b).  In  contrast  to  the  case  of 
the  transmission  coefficient,  the  reflection  does  not  admit  a  simple  formula  in  terms 
of  the  eigenvalues  of  the  system.  This  is  only  possible  when  the  Green’s  function  is 
separable,  namely,  when  the  distance  between  the  atoms  is  d  =  n\p/2. 


4-3.  Photonic  crystal  bandgap s 


The  band-gap  region  of  a  photonic  crystal  waveguide  (PCW)  is  a  very  appealing  scenario 
to  explore  coherent  atom-atom  interactions,  as  light  cannot  propagate,  and  atoms 
interact  with  each  other  through  evanescent  fields  30  .  For  a  photonic  crystal  waveguide 


of  lattice  constant  a  the  elements  of  matrix  g  are  well  approximated  by 

Qij  =  Jib  cos(nXi/a)  cos^TTXj /a)e~Kx^Xi~x^ , 


(22) 


where  the  cosine  terms  account  for  the  spatial  profile  of  the  Bloch  modes,  and  n~  1  is 
the  finite  range  of  interaction  due  to  the  evanescent  decay  of  the  guided  mode  field  in 
the  bandgap,  which  is  controlled  by  detuning  the  band-edge  frequency  from  the  atomic 
resonance.  It  should  be  noted  that  in  this  idealized  picture,  gtj  is  purely  real,  indicating 
the  absence  of  collective  emission  into  the  PCW.  This  is  naturally  expected,  due  to  the 
absence  of  guided  modes  at  the  atomic  frequency.  In  practice,  residual  decay  might  still 
exist  to  the  extent  that  the  mediating  photon  has  a  decay  channel.  This  could  be  either 
due  to  the  finite  length  of  the  PCW,  which  can  cause  the  photon  to  leak  out  the  ends 
and  is  suppressed  when  kx  L  g>>  1,  or  through  scattering  and  absorption  losses  of  the 
PCW.  Given  that  these  photonic  decay  processes  can  be  made  small,  for  conceptual 
simplicity  here  we  treat  the  idealized  case. 

For  a  chain  of  periodically  spaced  atoms  placed  in  even  antinodes  of  the  Bloch 
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modes,  the  dipole- projected  Green’s  function  matrix  reads 


0 


Ji 


1 

X 


ID 


x  x2 

1  X 


x^ 

xN~2 


(23) 


Vx^"1  x^-2  x^"3  •••  i  / 


where  we  have  defined  x  =  e~Kxd,  with  d  being  the  distance  between  nearest-neighbor 
atoms.  The  matrix  g  is  a  real  symmetric  Toeplitz  matrix  (or  bisymmetric  matrix). 
Neglecting  higher  order  contributions  besides  first- neighbor,  an  approximation  valid  for 
nxd  1,  g  becomes  a  tridiagonal  Toeplitz  matrix  whose  eigenvalues  and  eigenvectors 
are  [6l] : 

\  =  Jid,?  =  1  +  2e~Kxd  cos  ,  (24a) 

=  /vTT sin  (Jfi)  '  <24b) 

In  this  simple  tight  binding  model,  the  frequency  shifts  of  the  collective  atomic 
modes  are  distributed  around  Jid  with  a  frequency  spread  controlled  by  kx  (i.e. ,  for 
larger  kx ,  the  modes  are  closer  in  frequency).  However,  if  the  interaction  length  is  very 
large  compared  to  the  distance  between  the  atoms,  the  approximation  of  neglecting 
higher  order  neighbors  falls  apart,  and  the  eigenvalues  start  to  show  a  different  behavior. 
Eventually,  when  the  interaction  length  becomes  infinite  (or  much  larger  than  the  length 
of  the  atomic  cloud),  there  is  only  one  bright  mode,  of  eigenvalue  Ab  =  iVJiD.  This  is 
analogous  to  the  cavity  case,  where  the  interaction  range  is  also  infinite,  except  now  the 
eigenvalue  is  purely  real.  This  can  be  observed  in  Fig.  [4j  which  shows  how  the  collective 
frequency  shifts  coalesce  towards  Jid  for  large  nxd.  The  band-edge  of  a  photonic  crystal 
is  thus  a  cross-over  region  in  which  the  single  bright  mode  approximation  holds  and 
then  transitions  to  another  regime  where  it  breaks  down,  as  the  guided  mode  becomes 
evanescent  and  decays  substantially  within  the  length  of  the  PCW.  Importantly,  the 
bandgap  of  a  photonic  crystal  provides  a  tunable  interaction  range,  a  feature  which  is 
unique  to  this  kind  of  nanostructure,  and  makes  PCWs  remarkably  different  reservoirs 
from  either  cavities  or  unstructured  waveguides. 

In  the  following  section,  we  present  some  predictions  for  the  transmission  spectrum 
of  two  atoms  coupled  to  a  PCW  for  r1D  and  Jid  values  that  can  be  achieved 
experimentally  in  the  coming  years.  We  hope  that  the  foreseen  large  coherent  couplings 
between  the  atoms  combined  with  low  dissipation  through  the  guided  mode  help  to 
stimulate  a  new  generation  of  experiments  that  go  beyond  the  current  state  of  the  art. 


5.  Experimental  perspectives 


In  a  recent  experiment  36  ,  the  authors  have  observed  signatures  of  collective  atom- 
light  interactions  in  the  transmission  spectra  of  atoms  coupled  to  an  alligator  photonic 
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Figure  5.  (a)  Magnitude  of  the  ratio  between  the  coherent  and  dissipative  couplings 
through  the  guided  mode  of  an  alligator  PCW  |36j.  The  dashed  line  shows  the  ratio 
as  given  in  Fig.4  of  Ref.  36  ,  and  the  continuous  curve  represents  the  expected  ratio 


that  could  be  achieved  within  the  next  years  (see  text  for  more  details),  (b)  Evolution 
of  the  excited  state  population  of  atom  1  (blue  curve)  and  2  (orange  curve)  after  fully 
inverting  atom  1  at  the  initial  time.  The  resonance  frequency  of  the  atoms  lies  in  the 
bandgap  of  the  photonic  crystal,  with  the  atoms  placed  at  successive  even  antinodes 
(continuous  curve).  The  dashed  line  represents  the  non-interacting  scenario,  where  the 
off-diagonal  terms  of  g  are  zero.  The  spin  exchange  and  decay  rates  are  chosen  to  be 
Jid  =  —  3r0,  Tid  =  0.15ro,  and  T'  =  0.5ro.  The  lattice  constant  is  a  =  370  nm  and 
the  range  of  interaction  is  n~l  =  80a. 


crystal  waveguide.  They  have  recorded  these  spectra  for  various  frequencies  around 
the  band  edge  of  the  PCW,  exploring  different  physical  regimes.  Outside  the  bandgap, 
due  to  the  finite  size  of  the  PCW,  they  observe  the  formation  of  a  low-finesse  cavity 
mode  [as  shown  in  Fig.  3(a)  of  Ref.  |36),  at  a  frequency  zq].  At  resonance  with  this 
cavity  mode,  the  dissipative  single-atom  coupling  to  the  structure  is  Pid(pl)  —  1.5r0, 
as  obtained  from  steady-state  transmission  lineshape  measurements.  The  decay  rate 
into  leaky  modes  is  T'/To  —  1.1,  estimated  from  finite-difference  time-domain  (FDTD) 
numerical  calculations. 

After  tuning  the  spectral  features  of  the  PCW  so  that  the  resonance  frequency  of  the 
atoms  moves  into  the  bandgap,  they  observe  asymmetric  lineshapes,  revealing  significant 
coherent  coupling.  Specifically  at  z/bg  =  60  GHz  inside  the  bandgap,  the  spin  exchange 
and  decay  rates  are  Jid(^bg)/To  ~  —0.2  and  TiD^BGVro  —  0.01,  respectively.  Due  to 
the  evanescent  character  of  the  field  in  the  bandgap,  the  interaction  range  is  finite,  and 
at  z'bg  its  value  is  k~ 1  ~  80a,  being  a  =  370  nm  the  lattice  constant  of  the  alligator 
PCW.  While  this  experiment  constitutes  the  first  observation  of  more  than  one  emitter 
interacting  through  the  guided  modes  around  the  band  edge  of  a  PCW,  the  values  of  Jid 
and  T id  are  not  yet  good  enough  to  observe  direct  signatures  of  atom-atom  interactions 
such  as  time-dependent  spin  exchange.  Nevertheless,  we  expect  that  near-term  advances 
of  the  current  set  up  will  yield  dramatic  improvements  on  these  rates,  opening  the  door 
to  exploring  exciting  collective  atomic  phenomena. 

In  particular,  instead  of  using  an  alligator  PCW,  one  can  employ  a  slot  photonic 
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crystal  waveguide  (4,62  ,  i.e.  a  quasi-ID  waveguide  embedded  in  a  2D  photonic  crystal. 
This  structure  would  be  advantageous  due  to  several  reasons.  First  of  all,  it  inhibits 
atomic  emission  into  non-guided  modes  due  to  the  surrounding  2D  photonic  bandgap 
that  reduces  the  modes  into  which  the  atom  can  radiate.  Absent  inhomogeneous 
broadening,  early  simulations  demonstrate  that  it  is  possible  to  achieve  a  very  small 
non-guided  decay  rate,  i.e.  T'  ~  0.5r0.  Moreover,  one  can  engineer  flatter  bands, 
which  leads  to  an  increase  of  the  group  index  of  ng  ~  30  near  the  band-edge  (three 
times  larger  than  that  of  the  current  alligator),  according  to  FDTD  simulations.  Then, 
both  Jid  and  Y  m  would  experience  a  three- fold  increase.  Finally,  by  trapping  the 
atoms  at  the  center  of  the  nanostructure,  in  between  the  two  slots  and  not  above  as 
it  is  currently  done,  we  have  estimated  that  Jid  and  r  1D  would  be  five  times  larger. 
Summarizing,  we  project  FidOg  ) /To  —  22  at  the  first  cavity  resonance.  This  yields  the 
values  of  JidTbgI/To  —  —3  and  FiD(z/Bg)/Fo  —  0.15  for  a  detuning  from  the  band  edge 
i'bg  =  20  GHz,  where  the  range  of  interaction  is  nf1  ~  80a. 

Figure  (5](a)  compares  the  ratio  Jid /Y id  between  the  coherent  and  dissipative 
guided- mode  rates  for  the  current  alligator  PCW  (dashed  fine)  and  the  described  slot 
PCW  (continuous  line).  The  improved  ratio  for  the  later  structure  can  already  be 
observed  at  frequencies  just  beyond  the  band-edge,  and  becomes  |  Jid/Tid|  —  104  at 
a  detuning  of  0.5  THz  from  the  band-edge.  An  indisputable  signature  of  collective 
behavior  is  represented  in  Fig.  lib)  ,  which  shows  the  evolution  of  the  excited  state 
populations  of  two  atoms  placed  at  successive  even  antinodes  (continuous  curve),  after 
initially  inverting  one  of  them.  The  atoms  interact  through  the  guided  modes  of  the 
already  described  slot  PCW,  and  their  resonance  frequency  lies  inside  the  bandgap,  at 
the  frequency  for  which  the  interaction  range  is  Kf1  ~  80a.  The  dashed  lines  show  the 
expected  result  for  non  interacting  atoms,  where  the  off-diagonal  terms  of  0  are  zero,  a 
situation  that  occurs  when  the  atoms  are  separated  by  a  distance  d  nf 1  • 

To  summarize,  we  believe  that  there  is  a  bright  future  for  experiments  involving  not 
only  atoms,  but  also  superconducting  qubits  interacting  through  the  guided  modes  of  a 
microwave  photonic  crystal.  In  a  recent  experiment,  a  ratio  of  Tid/T'  =  50  has  already 
been  achieved  for  transmon  qubits  connected  to  a  ID  coplanar  microwave  transmission 
line  1 23).  Combined  with  the  exciting  recent  advances  in  microwave  photonic  crystal 


fabrication  3?  ,  we  expect  a  next  generation  of  experiments  where  many  qubits  interact 
with  each  other  in  a  mostly  coherent  manner. 


6.  Conclusion 

We  have  analyzed  the  optical  response  of  a  chain  of  atoms  placed  along  a  quasi-ID 
nanophotonic  structure  in  terms  of  the  classical  electromagnetic  Green’s  function.  This 
formalism  is  valid  in  the  presence  of  absorptive  and  dispersive  media. 

We  find  that  the  linear  response  of  the  atoms  can  be  understood  in  terms  of 
collective  atomic  eigenstates  of  the  Green’s  function  matrix  g(ay,  xf)  for  all  pairs  of 

atoms.  In  particular,  we  have  derived  a  closed  expression  for  the  transmission  spectra 
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that  only  depends  on  the  cooperative  frequency  shifts  and  decay  rates  of  these  modes. 
We  have  shown  that  the  transmission  coefficient  is  a  direct  probe  of  the  Green’s  function 
of  the  nanostructure,  enabling  us  to  determine  whether  the  atom-light  interactions  are 
fundamentally  dispersive  or  dissipative  in  character  as  well  as  to  quantify  the  degree  of 
cooperative  interaction.  We  have  gained  insight  into  the  interactions  between  atoms  and 
quasi-ID  cavities,  waveguides,  and  photonic  crystals,  structures  of  relevance  in  recent 
experiments,  as  well  as  provided  estimations  of  what  can  be  observed  in  the  near  future. 

The  Green’s  function  formalism  provides  a  natural  language  that  unifies 
nanophotonics  and  quantum  optics,  and  our  results  apply  not  only  to  atoms  36  ,  but 


to  many  other  quantum  emitters,  such  as  superconducting  qubits  |37  ,  NV  centers  63 


rare  earth  ions  [64|  or  quantum  dots  |4],  interacting  with  any  kind  of  quasi-ID  photonic 
structures  or  circuits. 


Appendix  A.  Transmission  and  reflection  coefficients  in  terms  of  Green’s 
functions 


We  begin  by  recalling  Eq.  (10), 


N 


E+(x)  =  Epx)  =  Epx)  -  y 


(gT(x)  -v£)  (q-Ej) 


£=i  (Aa  +  >A.  id)  +  i(F  + 

which  relates  the  field  along  any  point  of  the  structure  with  the  collective  atomic  modes. 
In  order  to  calculate  the  transmission  spectra,  we  need  an  expression  that  connects  the 
output  and  the  input  fields.  To  do  so,  let’s  consider  that  we  have  a  dipole  pieft  placed  to 
the  left  of  the  first  atom  of  the  chain,  at  position  a:ieft,  which  is  the  source  of  the  probe 
field  Ep  .  For  the  sake  of  simplicity,  pieft  is  polarized  along  y,  the  same  polarization 
of  the  guided  mode  field.  To  obtain  the  transmission  coefficient,  we  evaluate  the  field 
at  position  xright,  immediately  to  the  right  of  the  last  atom  of  the  chain.  When  the 
atoms  are  not  present,  the  probe  field  at  the  left  and  right  positions  of  the  quasi-ID 
nanostructure  are 

Ep  (fCleft)  /JjQUJp  GiD^ieft,  fCleft)  Pleft,  (A. la) 

Ep  (bright)  bO^'p  Gid (bright,  Teft)  Pleft-  (A. lb) 

Then,  the  transmission  for  the  system  without  the  atoms  is  simply 


*o(Aa)  — 


Ep  (bright)  _  Giu  (xright ,  ) 

ET(xieft)  G\Y) (^Cleft ,  ^left) 


(A.2) 


When  N  atoms  are  placed  in  the  vicinity  of  the  nanostructure,  the  field  at  position 

bright  IS 

1 


E  (xright)  Ep  (bright) 


=  t0(AA)  ~ 


y.  (g^right)  •  Vg)  (v[  '  g(^left)) 
g{x left)  -7’ left )  £=1  (Aa  +  ^  id)  +  i(T/  +  T^id)/2  P 

1 


^(Xieft,  3?  left)  t=1  (Aa  +  Id)  +  l(T '  +  1d)/2 
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where  we  have  employed  that  the  probe  field  at  atom  Xj  can  be  related  to  E+  (xieft) 
as  E+(Xj)  =  [IqOJ2  Gid  (Xj ,  -Tieft)  Pleft  =  Ep  •  Then>  the  normalized 


transmission  coefficient  is 

*(Aa)/*o(Aa)  =  1  - 


1 


^  (gT(^right)  '  V$)  (v[  •  g(xleft))  /  A  ^ 

7a  ;  T  TTTrFv- wo’ 


9  (bright  5  %  left)  £_  i  (^A  "t"  «^£,1d)  “1“  f(T  ~L  r^,lD)/2 

as  shown  in  the  main  text.  Let’s  now  calculate  the  reflection  coefficient.  Without  the 
atoms,  the  field  at  x\eft  is  E+(x\ett)  =  [l+r0(AA)]Lj"(rieft).  When  the  atoms  are  present, 
the  field  reads 


E+(xieh)  =  [1  +r0(AA)]£p  (iCleft) 


1 


g(xieft,  ^ieft)  ^=1  (Aa  +  Aid)  +  i(rf  +  rCiiD)/2 

Following  similar  steps  as  those  above,  we  find 

N 


r{ Aa)  =  r0(AA)  - 


1 


£  (gT(^left)  •  Vg)  (v[  '  g(^left))  (A  4) 


g(x left)  left )  £=1  (Aa  +  J^,id)  +  i(r'  +  r€;1D)/2’ 


the  equation  in  the  main  text. 


Appendix  B.  Derivation  of  Equation  (12)  for  the  transmission 


We  can  exploit  some  properties  of  ID  systems  to  arrive  to  the  closed  expression  for  the 


transmission  shown  in  Eq.  (12),  which  only  depends  of  the  decay  rates  and  frequency 


shifts  of  the  modes,  not  on  their  spatial  structure  (i.e.  the  eigenfunctions).  We  first 
show  how  to  derive  the  ID  Green’s  function  wave  equation,  and  how  the  solution  is 
related  to  the  full  quasi-ID  solution.  We  start  with  the  3D  Green’s  function  Gid  for  the 
guided  mode,  which  follows  Eq.  Q.  We  assume  that  the  guided  modes  are  transverse 
waves  that  travel  in  the  ±x  direction  and  are  polarized  along  y,  and  that  the  field 
is  approximately  uniform  in  the  transverse  directions.  From  3D,  one  can  in  principle 
construct  the  guided  modes  and  their  dispersion  relations  u(k),  from  which  one  can 
identify  an  effective  dielectric  constant  eeff(x,u>)  which  produces  the  same  behavior  (at 
least  within  some  bandwidth).  The  final  answer  that  we  are  trying  to  achieve  does  not 
depend  on  explicit  construction  of  eeff(x,(u).  The  result  is  a  Helmoltz  equation  for  the 
Green’s  function  that  reads 


d2u2 

-  +  -eesM 


Gid(x,x'  ,u)  =  —8(x  —  x1), 


(B.l) 


where  Gid  =  AG  id,  being  A  the  effective  mode  area.  The  solution  for  this  second  order 
linear  ordinary  differential  equation  can  be  expressed  as  the  sum  of  the  two  homogeneous 
solutions.  The  Green’s  function  can  then  be  written  in  terms  of  the  auxiliary  fields 
4>l,r(x),  which  are  solutions  of  the  homogeneous  equation,  as 


G1D(x,x')  = 


Q(x'  -  x)4>l(x')4>r(x)  +  Q(x  -  x')4>R{x’)f)L{x ) 
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where  W  is  the  Wronskian,  which  does  not  depend  on  the  position,  and  is  given  by 

(B.3) 


Tjr/  A  Mx')  dfR(x,)1  , 

W  =  <f>R(x  )— — - —  (Pl{x  ). 


dx'  dx‘ 

We  can  then  recover  the  full  Green’s  function  between  atom  i  and  atom  j  as 
Ginixi,  Xj,uj)  =  ^-G1D(xi,Xj,uj)  =  [Q(xj  -  Xi)<j>L(xj)<j>R{xi)  +  <d(xi  -  xi)0L(^)0R(xi)] , 

(B.4) 

where  <f>R =  4>r,l/'/AW  .  Then,  the  dipole-projected  Green’s  function  is 

9ij  =  ~  Xi)sji  -  (-)(.r,  -  Xj)sij, 


where  Sij  =  (pL(xi)(pR(xj),  with  cpLti  =  ^JJl0uJpd?/h  (f>L(xi)  and  (pRJ  = 


PoiOpd2 /h  (pR(x:j).  It  is  convenient  to  define  the  rank-one  matrix  s  =  tpi  0  where 
^P{r,l}  =  (<P{r,l}(x i), (P{r,l}(xn))  is  a  vector  of  N  components.  Let’s  now  proceed  to 
demonstrate  Eq.  (12).  In  terms  of  the  eigenfunctions  of  0,  the  transmission  is 

(gT  (fright)  '  V$)  (vj  ■  g(xieft)) 


t(AA)/t0(AA)  =  1  - 
=  1  - 


N 

E 


9 (bright ?  % left)  £— \  (^A  T  IE))  T  i(T  ~\~  T^5Xd)/2 
9fe8h.,x,eft)  (gT(Xrlgh,)  '  M~1  '  e(XkG  > 


where  M  is  given  in  Eq.  (|8j).  Since  g  oc  Gid,  and  using  the  expression  for  the  Green’s 
function  in  terms  of  the  right-going  and  left-going  field  solutions  [Eq.  (B.4)],  we  find 

1  1 


t(AA)/t0(AA)  =  l-(pR- 


Aa  +  iP/2  +  g 


•<Pl  =  1  ~  v 


1  +  0 


U, 


where  we  have  defined  v  =  ipR/y  Aa  +  iP/2,  u  =  <£l/ \J Aa  +  iP/2,  and  0  =  0/(Aa  + 
iT;/2).  By  the  matrix  determinant  lemma  (54] ,  we  know  that  for  a  invertible  matrix  A 
and  a  pair  of  vectors  u,  v,  we  can  write  det(A  +  u  ®  vT)  =  det(A)  (l  +  vT  •  A-1  •  u) . 
Choosing  A  =  —  (1  +  0),  we  find 


f(AA)/f0(AA)  = 


det(l  +  0  —  u  ®  vT)  det((AA  +  iT'/2)l  +  0  —  5) 


det(l  +  0)  det((AA  +  iP/2)!  +  0) 

Since  (Aa  +  ir'/2)l  +  0  —  5  is  a  triangular  matrix  with  (Aa  +  iP/2)  in  the  diagonal 
entries,  and  the  determinant  of  a  triangular  matrix  is  the  product  of  the  diagonal  entries, 
we  find  det((AA  +  ir'/2)l  +  0  —  s)  =  (Aa  +  iT'/2)N,  which  yields 


t(AA)/t0(AA)  = 


(Aa  +  ir'/2) 


N 


N 

n 


Aa  +  iP/2 


det((AA  +  ir,/2)3L  +  g)  ,=1  (Aa  +  id)  +  i(r'  +  T^id)/2 


as  the  determinant  of  a  matrix  is  the  product  of  its  eigenvalues.  The  above  expression  is 

precisely  Eq.  (ig.  To  the  best  of  our  knowledge,  it  is  not  possible  to  obtain  a  simplified 

expression  for  the  reflection  coefficient.  .  ... 
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We  present  a  platform  for  the  simulation  of  quantum  magnetism 
with  full  control  of  interactions  between  pairs  of  spins  at  arbitrary 
distances  in  ID  and  2D  lattices.  In  our  scheme,  two  internal  atomic 
states  represent  a  pseudospin  for  atoms  trapped  within  a  photonic 
crystal  waveguide  (PCW).  With  the  atomic  transition  frequency 
aligned  inside  a  band  gap  of  the  PCW,  virtual  photons  mediate 
coherent  spin-spin  interactions  between  lattice  sites.  To  obtain  full 
control  of  interaction  coefficients  at  arbitrary  atom-atom  separa¬ 
tions,  ground-state  energy  shifts  are  introduced  as  a  function  of 
distance  across  the  PCW.  In  conjunction  with  auxiliary  pump  fields, 
spin-exchange  versus  atom-atom  separation  can  be  engineered 
with  arbitrary  magnitude  and  phase,  and  arranged  to  introduce 
nontrivial  Berry  phases  in  the  spin  lattice,  thus  opening  new  ave¬ 
nues  for  realizing  topological  spin  models.  We  illustrate  the  broad 
applicability  of  our  scheme  by  explicit  construction  for  several  well- 
known  spin  models. 

nanophotonics  |  quantum  matter  |  cold  atoms  |  quantum  many-body  | 
quantum  spin 

Quantum  simulation  has  become  an  important  theme  for  re¬ 
search  in  contemporary  physics  (1).  A  quantum  simulator 
consists  of  quantum  particles  (e.g.,  neutral  atoms)  that  interact  by 
way  of  a  variety  of  processes,  such  as  atomic  collisions.  Such  pro¬ 
cesses  typically  lead  to  short-range,  nearest-neighbor  interactions 
(2-6).  Alternative  approaches  for  quantum  simulation  use  dipolar 
quantum  gases  (7,  8),  polar  molecules  (9-11),  and  Rydberg  atoms 
(12-15),  leading  to  interactions  that  typically  scale  as  1/r3,  where  r  is 
the  interparticle  separation.  For  trapped  ion  quantum  simulators 
(16-20),  tunability  in  a  power  law  scaling  of  r~}]  with  0<rj<3  can  in 
principle  be  achieved.  Beyond  simple  power  law  scaling,  it  is  also 
possible  to  engineer  arbitrary  long-range  interactions  mediated  by 
the  collective  phonon  modes,  which  can  be  achieved  by  independent 
Raman  addressing  on  individual  ions  (21). 

Using  photons  to  mediate  controllable  long-range  interactions 
between  isolated  quantum  systems  presents  yet  another  approach 
for  assembling  quantum  simulators  (22).  Recent  successful  ap¬ 
proaches  include  coupling  ultracold  atoms  to  a  driven  photonic 
mode  in  a  conventional  mirror  cavity,  thereby  creating  quantum 
many-body  models  (using  atomic  external  degrees  of  freedom) 
with  cavity-field-mediated  infinite-range  interactions  (23). 
Finite-range  and  spatially  disordered  interactions  can  be  realized 
by  using  multimode  cavities  (24).  Recent  demonstrations  on 
coupling  cold  atoms  to  guided  mode  photons  in  photonic  crystal 
waveguides  (25,  26)  and  cavities  (27,  28)  present  promising  ave¬ 
nues  (using  atomic  internal  degrees  of  freedom)  due  to  unprec¬ 
edented  strong  single  atom-photon  coupling  rate  and  scalability. 
Related  efforts  also  exists  for  coupling  solid-state  quantum  emit¬ 
ters,  such  as  quantum  dots  (29,  30)  and  diamond  nitrogen-vacancy 
centers  (31,  32),  to  photonic  crystals.  Scaling  to  a  many-body 
quantum  simulator  based  on  solid-state  systems,  however,  still 
remains  elusive.  Successful  implementations  can  be  found  in 
the  microwave  domain,  where  superconducting  qubits  behave  as 
artificial  atoms  strongly  coupled  to  microwave  photons  propa¬ 
gating  in  a  network  formed  by  superconducting  resonators  and 
transmission  lines  (33-35). 


Here,  we  propose  and  analyze  a  physical  platform  for  simulating 
long-range  quantum  magnetism  in  which  full  control  is  achieved  for 
the  spin-exchange  coefficient  between  a  pair  of  spins  at  arbitrary 
distances  in  ID  and  2D  lattices.  The  enabling  platform,  as  described 
in  refs.  36  and  37,  is  trapped  atoms  within  photonic  crystal  wave¬ 
guides  (PCWs),  with  atom-atom  interactions  mediated  by  photons 
of  the  guided  modes  (GMs)  in  the  PCWs.  As  illustrated  in  Fig.  1 A 
and  B,  single  atoms  are  localized  within  unit  cells  of  the  PCWs  in  ID 
and  2D  periodic  dielectric  structures.  At  each  site,  two  internal 
atomic  states  are  treated  as  pseudospin  states,  with  spin- 1/2  con¬ 
sidered  here  for  definiteness  (e.g.,  states  \g)  and  | s)  in  Fig.  1C). 

Our  scheme  uses  strong,  and  coherent  atom-photon  interactions 
inside  a  photonic  band  gap  (36-40),  and  long-range  transport 
property  of  GM  photons  for  the  exploration  of  a  large  class  of 
quantum  magnetism.  This  is  contrary  to  conventional  hybrid 
schemes  based  on,  for  example,  arrays  of  high  finesse  cavities 
(41-44)  in  which  the  pseudospin  acquires  only  the  nearest  (or  at 
most  the  next-nearest)  neighbor  interactions  due  to  strong  expo¬ 
nential  suppression  of  photonic  wave  packet  beyond  single  cavities. 

In  its  original  form  (36-40),  the  localization  of  pseudospin  is 
effectively  controlled  by  single-atom  defect  cavities  (36).  The  cavity 
mode  function  can  be  adjusted  to  extend  over  long  distances  within 
the  PCWs,  thereby  permitting  long-range  spin  exchange  interac¬ 
tions.  The  interaction  can  also  be  tuned  dynamically,  via  external 
addressing  beams,  to  induce  complex  long-range  spin  transport, 
which  we  describe  in  the  following  (36,  37). 

To  engineer  tunable,  long-range  spin  Hamiltonians,  we  use  an 
atomic  A  scheme  and  two-photon  Raman  transitions,  where  an  atom 
flips  its  spin  state  by  scattering  one  photon  from  an  external  pump 
field  into  the  GMs  of  a  PCW.  The  GM  photon  then  propagates 
within  the  waveguide,  inducing  spin  flip  in  an  atom  located  at  a 
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Fig.  1.  Photon-mediated  atom-atom  interactions  in  (A)  ID  and  (B)  2D  PGA/s. 
(C)  Atomic-level  scheme:  atomic  dipole  |s)  <-»•  |e)  is  coupled  to  an  external  pump, 
|g)  <-►  |e)  coupled  to  a  GM  photon,  and  r*,  the  excited  state  decay  rate  to  free  space 
and  leaky  modest  (D)  Simplified  band  structure  ®(k)  near  the  band  edge  k  =  kc  and 
&>(kc)  =a>c-  Atomic  transition  frequency  (oeg  =  (oe-(og  lies  within  the  band  gap. 


distant  site  via  the  reverse  two-photon  Raman  process.  When  we 
align  the  atomic  resonant  frequency  inside  the  photonic  band  gap,  as 
depicted  in  Fig.  ID,  only  virtual  photons  can  mediate  this  remote 
spin  exchange  and  the  GM  dynamics  are  fully  coherent,  effectively 
creating  a  spin  Hamiltonian  with  long-range  interactions.  As  dis¬ 
cussed  in  refs.  36  and  37,  the  overall  strength  and  length  scale  of  the 
spin-exchange  coefficients  can  be  tuned  by  an  external  pump  field, 
albeit  within  the  constraints  set  by  a  functional  form  that  depends  on 
the  dimensionality  and  the  photonic  band  structure.  These  con¬ 
straints  may  limit  our  ability  to  explore  novel  quantum  phases  and 
nonequilibrium  dynamics  in  various  spin  models,  because  many  ef¬ 
fects  display  strong  dependencies  on  the  functional  form  of  long- 
range  interactions  (45-50).  It  is  therefore  highly  desirable  to  obtain 
full  control  of  interactions  without  the  need  to  investigate  over  a 
wide  range  of  PCW  designs  with  different  photonic  band  structures. 

To  fully  control  spin-exchange  coefficients  at  arbitrary  separations, 
here  we  adopt  a  Raman-addressing  scheme  similarly  discussed  for 
cold  atoms  and  trapped  ions  (51-55).  We  introduce  atomic  ground- 
state  energy  shifts  as  a  function  of  distance  across  the  PCW.  Due  to 
conservation  of  energy,  these  shifts  suppress  reverse  two-photon 
Raman  processes  in  the  original  scheme  (36,  37),  forbidding  spin 
exchange  within  the  entire  PCW.  However,  we  can  selectively  acti¬ 
vate  certain  spin-exchange  interactions  J(rm/l )  between  atom  pairs 
(m,n)  separated  by  rm?n,  by  applying  an  auxiliary  sideband  whose 
frequency  matches  that  of  the  original  pump  plus  the  ground-state 
energy  shift  between  the  atom  pairs.  This  allows  us  to  build  a  pre¬ 
scribed  spin  Hamiltonian  with  interaction  terms  “one  by  one.”  Note 
that  each  sideband  in  a  Raman-addressing  beam  can  be  easily  in¬ 
troduced,  for  example,  by  an  electro-optical  modulator.  By  in¬ 
troducing  multiple  sidebands  and  by  controlling  their  frequencies, 
amplitudes,  and  relative  phases,  we  can  engineer  spin  Hamiltonians 
with  arbitrary,  complex  interaction  coefficients  J(rm/l).  Depending 
on  the  dimensionality  and  the  type  of  spin  Hamiltonians,  our  scheme 
requires  only  one  or  a  few  Raman  beams  to  generate  the  desired 
interactions.  Furthermore,  by  properly  choosing  the  propagation 
phases  of  the  Raman  beams,  we  can  imprint  geometric  phases  in  the 
spin  system,  thus  providing  unique  opportunities  for  realizing  topo¬ 
logical  spin  models. 

We  substantiate  the  broad  applicability  of  our  methods  by 
explicit  elaboration  of  the  set  of  pump  fields  required  to  realize 
well-known  spin  Hamiltonians.  For  ID  spin  chains,  we  consider 
the  implementation  of  the  Haldane-Shastry  model  (56,  57).  For 
2D  spin  lattices,  we  elaborate  the  configurations  for  realizing 
topological  flat  bands  (58,  59)  in  Haldane’s  spin  model  (56),  as 
well  as  a  “checkerboard”  chiral-flux  lattice  (58,  59).  We  also 
consider  a  2D  XXZ  spin  Hamiltonian  with  J(rm,n)  <x  l/r^n  and 


rj  =  1,2,3  (60).  In  addition,  we  report  numerical  results  on  the  rj 
dependence  of  its  magnetization  diagram. 

Controlling  Spin-Spin  Interaction  Through  Multifrequency 
Driving 

In  the  following,  we  discuss  how  to  achieve  full  control  of  interac¬ 
tions  by  multifrequency  pump  fields.  We  assume  (i)  N  atoms  trapped 
in  either  a  ID  or  2D  PCW,  as  depicted  in  Fig.  1  A  and  B,  with  a 
spatially  dependent  ground-state  energy  shift  oog.  For  simplicity,  we 
assume  one  atom  per  unit  cell  of  the  PCW,  although  this  assumption 
can  be  relaxed  afterward;  (ii)  the  structure  is  engineered  (22-28) 
such  that  the  GM  polarization  is  coupled  to  the  atomic  dipole, 

| g)  ++  \e),  as  shown  in  Fig.  1C,  and,  under  rotating  wave  approxi¬ 
mation,  is  described  by  the  following  Hamiltonian  (using  h  =  1): 

Hi  m  =  ^2gk(rn)ako,‘g  +  h.c.,  [1] 

M 

where  gk(rn)  =gkelk  Tn  is  the  single-photon  coupling  constant  at 
site  location  rn,  with  n  being  the  site  index;  ak,  the  GM  field 
operator;  and  (Aab  =  \a)n(b\,  the  atomic  operators  with  a,b  being 
one  of  the  g,s,e  states.  Moreover,  as  in  refs.  36  and  37,  we 
assume  (Hi)  there  is  another  hyperfine  level  | s),  addressed  by  a 
Raman  field  with  coupling  strength  Q  as  follows: 

Hd  W=E(£y^/w+hx)>  [2] 

where  col  is  the  main  driving  frequency.  The  Raman  field  Q(t) 
contains  mp  frequency  components  that  are  introduced  to 
achieve  full  control  of  the  final  effective  spin  Hamiltonian.  Full 
dependence  of  Q(t)  can  be  written  as  follows: 

Hip- 1 

0(0  =  X!  [3] 

a— 0 

where  doa  are  the  detunings  of  the  sidebands  from  the  main  frequency 
coL  such  that  do o  =  0,  and  Qa,  the  complex  amplitudes. 

We  can  adiabatically  eliminate  the  excited  states  \e)  and  the 
photonic  GMs  under  the  condition  that  (iv)  max{|Q|,  \doa- do p\}  <C 
|  A |  =  \ooe  -  &>l|-  This  condition  guarantees  that,  first,  the  excited 
state  is  only  virtually  populated,  and  that,  second,  the  time  de¬ 
pendence  induced  by  the  sideband  driving  is  approximately  constant 
over  the  timescale  A-1.  As  discussed  in  refs.  36  and  37,  if  qol  -  aog 
lies  in  the  photonic  band  gap,  photon-mediated  interactions  by  GMs 
are  purely  coherent^  Under  the  Bom-Markov  approximation,  we 
then  arrive  at  an  effective  AT  Hamiltonian  (SI  Appendix  A:  Complete 
Derivation  of  Final  Time-Dependent  Hamiltonian)'. 

N  rap-1 

Hxy(t)=  Y,  [4] 

m,  n^m  a,  /i=0 

where  we  have  defined  Xa  =  Qa/(2A);  cog^n=cog( rn)  is  the  site- 
dependent  ground-state  energy  shift,  and  J(rm/l)  is  the  atom- 
GM  photon  coupling  strength  (36,  37)  that  typically  depends 
on  atomic  separation  rm  „  =rm-rn. 

We  focus  on  “sideband  engineering”  and  treat  J(rm/l)  as  ap¬ 
proximately  constant  over  atomic  separations  considered.* *  This 
is  valid  as  long  as  the  farthest  atomic  separation  with  nonzero 
engineered  interaction  is  much  smaller  than  the  decay  length 


+To  simplify  the  discussion,  in  this  paper,  we  neglect  decoherence  effects  caused  by 
atomic  emission  into  free  space  and  leaky  modes  as  well  as  photon  loss  due  to  imper¬ 
fections  in  the  PCW.  These  effects  were  both  carefully  discussed  in  refs.  36  and  37, 
suggesting  the  number  of  spin-exchange  cycles  in  the  presence  of  decoherence  can  re¬ 
alistically  reach  yf’»35~100  using  ultra-high  Q  PCWs. 

*One  may  also  replace  a  PCW  with  a  single-mode  nanophotonic  cavity,  operating  in  the 
strong  dispersive  regime  (61,  62),  to  achieve  constant  GM  coupling  J  independent  of 
|rm#n|.  Realistic  nanophotonic  cavity  implementations  will  be  considered  elsewhere. 
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[8] 


\e) 


resonant  Raman-scattering  processes.  Spin  exchanges  (/\)  \sn,gm)  -*■  \gn,sm) 
and  ( B )  \gn,  sm)  ->  |sn,  gm)  are  allowed  only  when  the  condition  - cogin  = 
cbp-cba  is  satisfied.  Qa/ A  and  Q^/A  control  the  exchange  rate. 


rap-1 

Jm,n  —  ^  ^  XaXpJ8n—m,fi—ce 
a,  (3=0 

It  can  be  shown  that  the  set  of  equations  defined  by  Eq.  8  has  at  least 
one  solution  for  any  arbitrary  choice  of  Jms,  that  is,  by  choosing 
£20  >  o  and  Jm,n  «  (XoX*_m  More  solutions  can 

be  found  by  directly  solving  the  set  of  nonlinear  equations  Eq.  8. 

It  is  important  to  highlight  that  multifrequency  driving  also 
enables  the  possibility  to  engineer  geometrical  phases  and, 
therefore,  topological  spin  models.  If  the  pump  field  propagation 
is  not  perfectly  transverse,  that  is,  kL  •  rm^  ±  0  (kL  being  the 
wave  vector  of  the  Raman  field),  the  effective  Hamiltonian  Eq.  7 
acquires  spatial-dependent,  complex  spin-exchange  coefficients 
via  the  phase  of  XaX ^  in  Eq.  8;  see  later  discussions. 

Beyond  an  ideal  setting,  we  now  stress  a  few  potential  error 
sources.  First,  for  practical  situations,  the  gradient  per  site  8  will 
be  a  limited  resource,  making  Eq.  7  not  an  ideal  approximation. 
Careful  Floquet  analysis  on  time-dependent  Hamiltonian  in  Eqs. 
5  and  6  is  required,  to  be  discussed  later.  Second,  there  is  an 
additional  Stark  shift  on  state  | s)  due  to  the  Raman  fields: 


scale  %=  of  the  coupling  strength  J(rm,n).  Here,  A  is  the 

band  curvature  (Fig.  ID),  Ac  =max{coc  -  (&>l  -<%„)}  is  the 
maximal  detuning  of  the  band  edge  to  the  frequency  of  coupled 
virtual  photons  that  mediate  interactions  (Fig.  1C),  and  we  have 
assumed  that  the  variation  of  ground-state  energies  cog,n  are  small 
compared  with  Ac.  Exact  functional  form  of  can  be  found 

in  refs.  36  and  37,  and  in  SI  Appendix  A:  Complete  Derivation  of 
Final  Time-Dependent  Hamiltonian. 

The  time  dependence  in  Eq.  4  can  be  further  engineered  and 
simplified.  We  note  that  the  interaction  between  two  atoms  n  and  m 
will  be  highly  dependent  on  the  resonant  condition  oog,m  -  cog^n  = 
mp  -  ooa,  provided  the  ground-state  energy  difference  \cog,n  -  cog/n  |  is 
much  larger  than  the  characteristic  timescale  of  interactions 
\XaX*J\.  The  intuitive  picture  is  depicted  in  Fig.  24:  the  atom 
n  scatters  from  sideband  a  a  photon  with  energy  qjl  +  ma  -  oog^n  into 
the  GMs.  When  this  GM  photon  propagates  to  the  atom  m,  it  will 
only  be  rescattered  into  a  sideband  f  that  satisfies  co^-\-coa- 
cogn  =(]0h  +  Q)p-  Q)gmj  whereas  the  rest  of  the  sidebands  remain 
off-resonant.  Fig.  2 B  depicts  a  reversed  process. 

For  concreteness,  we  discuss  a  ID  case  where  we  assume  (v)  a 
linear  gradient  in  the  ground-state  energy  cog^n  =  n8,  with  8  being 
the  energy  difference  between  adjacent  sites.  The  sidebands  will 
be  chosen  accordingly  such  that  cba  =  a8,  with  aG  Z. 

Summing  up,  with  all  these  assumptions  (i-v),  the  resulting 
effective  Hamiltonian  Eq.  4  can  finally  be  rewritten  as  follows: 

Hxy  (t )  =  Hxy,  p  eipdt,  [5] 

p 

where  Hxy,p  is  the  contribution  that  oscillates  with  frequency  p8. 
Written  explicitly, 


N  rap-1 

Hxy,p  =  y!  XaXp8n (/^s (Asg .  [6] 

m,n^m  a, (3=0 


In  an  ideal  situation,  the  gradient  per  site  satisfies  8  »  \xax;j\ 
such  that  the  contributions  from  Hxy,p  V/?/0  can  be  neglected. 
Under  these  assumptions,  we  arrive  at  an  effective  time-indepen¬ 
dent  Hamiltonian: 


N 

HXY(t)  ~HxY,Q  =  [7] 

m,n^m 


where  couplings  Jm ^  can  be  tuned  by  adjusting  the  amplitudes 
and  phases  of  the  sidebands  Xa  as  they  are  given  by  the  following: 


8cos(t)  =  - 


rap-l 

E 


a= 0 


|£4|2 

4A 


rap-1 

E* 


a>(3 


i(ma-Q)B)t 

2A 


[9] 


where  $R[.]  indicates  real  part.  We  note  that  the  time-independent 
contribution  in  Eq.  9  can  be  absorbed  into  the  energy  of  cos  without 
significant  contribution  to  the  dynamics,  whereas  the  time-depen¬ 
dent  terms  may  be  averaged  out  over  the  atomic  timescales  that  we 
are  interested  in.  We  will  present  strategies  for  optimizing  the 
choice  of  8,  and  minimizing  detrimental  effects  due  to  undesired 
time-dependent  terms  in  Eqs.  5  and  9  in  later  discussions. 

Independent  Control  of  XX  and  YY  Interactions.  So  far,  we  can  fully 
engineer  an  XY  Hamiltonian  with  equal  weight  between  XX  and 
YY  terms  by  defining  the  Pauli  operators  (px,  oy,  of)  —  (osg  + 
<7gs,  i(osg  -ogs),  ogg  -  oss).  We  now  show  flexible  control  of  AW  and 
YY  interactions  with  slight  modifications  in  the  atomic  level  structure 
and  the  Raman-addressing  scheme.  In  particular,  we  use  a  butterfly¬ 
like  level  structure  where  there  are  two  transitions,  \g)  \e)  and 

| s)  \e),  coupled  to  the  same  GM,  as  depicted  in  Fig.  3.  We  will  use 
two  multifrequency  Raman  pump  fields,  Qg(t)  and  Ds(t),  to  induce 
| g)  ++  \e)  <-»  | s)  and  | s)  \e)  ++  |g)  two-photon  Raman  transitions, 
respectively. 

For  example,  to  control  XX  or  YY  interactions,  we  require  that 
the  two  pump  fields  induce  spin  flips  with  equal  amplitude,  that 
is,  ogs  ±  osg.  This  is  possible  if  we  choose  the  main  frequencies  of 
the  pumps  (copg  and  cops)  such  that  cops  =  a>Lg  +  2cog,  and  match 
their  amplitudes  such  that  |£2g;Cr|/Ag  =  |£\a|/A5,  where  A5  = 
coe  ~  ^.g  —  We  ~  (^Ljg  T  (fig)->  &nd  |  A5jg| 

Adiabatically  eliminating  the  excited  states  as  well  as  the  GMs, 
we  arrive  at  the  following  Hamiltonian: 

Hxx,YY,0  =  £  [Jm,n  (og  +  ^4)  (<,  +  +  h.C.  ], 

m,n>m 

[10] 

where  (p^  is  the  relative  phase  between  the  pumps  fields  £2^.  Assum¬ 
ing  the  laser  beams  that  generate  the  Raman  fields  are  copropagating 
or  are  both  illuminating  the  atoms  transversely,  that  is,  kL  •  r m,n  =  0, 
we  can  generate  either  X  or  Y components,  (ofg  ±o^),  by  setting  the 
phase  (pgs  =  0  or  n;  more  exotic  combinations  are  available  with 
generic  choice  of  (pgs.  Moreover,  if  the  laser  beams  are  not  cop¬ 
ropagating,  they  create  spatially  dependent  phases  (pgs  m .  This  can 
create  site-dependent  AX,  YY,  or  XY  terms. 
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Fig.  3.  Atomic  "butterfly"  level  structure.  Two  pump  fields  Qs  and  Qg, 
tuned  to  couple  to  the  same  GM  photon,  are  introduced  to  control  XX  and 
YY  interactions  independently. 


Independent  Control  of  ZZ  Interactions.  An  independently  con¬ 
trolled  ZZ  Hamiltonian,  in  combination  with  arbitrary  XY  terms, 
would  allow  us  to  engineer  SU(2)-invariant  spin  models  as  well 
as  a  large  class  of  XXZ  models,  that  is,  the  following: 

HXXZ=HXY+HZZ=  jr  [(2C,„«  +  h.c.)+4X^]- 

m,n>m 

HU 


In  refs.  36  and  37,  it  was  shown  that  ZZ  interaction  can  be  created 
by  adding  an  extra  pump  field  to  the  | g)  ++  \e)  transition  in  Fig.  1C. 
However,  as  ZZ  terms  in  this  scheme  (36,  37)  do  not  involve  flipping 
atomic  states,  it  is  not  directly  applicable  to  our  multifrequency  pump 
method.  Nonetheless,  because  we  can  generate  XX  and  YY  inter¬ 
actions  independently,  a  straightforward  scheme  to  engineer  Hzz  is 
to  use  single  qubit  rotations  to  rotate  the  spin  coordinates  X++Z  or 
y^Z,  followed  by  stroboscopic  evolutions  (63)  to  engineer  the  full- 
spin  Hamiltonian.  Spin-rotation  can  be  realized,  for  example,  with 
a  collective  microwave  driving  Hmw  =  J]((£2mw/ 2)o^  +  h.c.),  in 

n 

which  a  ^/2-microwave  pulse  rotates  the  basis  |s)„j 

Thus,  an  Hxxz  Hamiltonian  can  be  simulated  using  the  fol¬ 
lowing  stroboscopic  evolution:  {Hxy,Hzz,Hxy,Hzz,  •  •  •  }  in  Nt 
steps  as  schematically  depicted  in  Fig.  4.  As  shown  in  SI  Appendix 
A:  Complete  Derivation  of  Final  Time-Dependent  Hamiltonian , 
the  error  accumulated  in  these  Nt  steps  can  be  bounded  by  the 
following: 


E2< 


N(RJt)2 


N, 


[12] 


where  J  =  max[/m  n]  is  the  largest  energy  scale  of  the  Hamiltonian  we 
want  to  simulate,  and  R  is  the  approximate  number  of  atoms  coupled 
through  the  interaction.  For  example,  if  Jm,n  is  a  nearest-neighbor 
interaction,  R  =  1.  If  Jm,n  oA/\m  -n\\  then  R  cx  Yln= l  V \n which 
typically  grows  much  slower  than  N.  Because  E2ocl  /Nt,  the  Trotter 
error  in  Nt  steps  can  in  principle  be  decreased  to  a  given  accuracy  s 
by  using  enough  steps,  that  is,  Nt  >  ( N(RJt)2)/s . 

More  complicated  stroboscopic  evolutions  may  lead  to  a  more 
favorable  error  scaling  (64-66),  although  in  real  experiments 
there  will  be  a  trade-off  between  minimizing  the  Trotter  error 
and  the  fidelity  of  the  individual  operations  to  achieve  Hxy  and 
Hzz.  As  this  will  depend  on  the  particular  experimental  setup,  we 
will  leave  such  analysis  out  of  current  discussions.  For  illustra¬ 
tion,  we  will  only  consider  the  simplest  kind  of  stroboscopic 
evolution  that  we  depicted  in  Fig.  4. 


Engineering  Spin  Hamiltonians  for  ID  Systems:  The 
Haldane-Shastry  S  =  1  /2  Spin  Chain 

In  the  first  example,  we  engineer  a  Haldane-Shastry  spin 
Hamiltonian  in  one  dimension  (56,  57): 


N-lN-m 

Hhs = z  z jn  [2  ( o^go^+n  +  h.c.)  +  ,  [13] 

m— 1  n- 1 

where  Jn  =Jo/sm2(nn/N),  Jq=Jji2  /TV2,  and  N  is  the  number  of 
spins.  The  interaction  strength  decays  slowly  with  approximately 
a  1/r2  dependence  while  satisfying  a  periodic  boundary  condi¬ 
tion.  Such  a  spin  Hamiltonian  is  difficult  to  realize  in  most  phys¬ 
ical  setups  that  interact,  for  example,  via  dipolar  interactions. 

We  can  engineer  the  periodic  boundary  condition  and  the 
long-range  interaction  Jn  directly  using  a  linear  array  of  trapped 
atoms  coupled  to  a  PCW.  To  achieve  this,  we  induce  atomic 
ground-state  energy  shift  mS  according  to  the  spin  index  m ,  and 
then  uniformly  illuminate  the  trapped  atoms  with  an  external 
pump  consisting  of  N  frequency  components  cba  =  ad ,  each  with 
an  amplitude  denoted  by  Qa  and  a  =  0,1,  . . .  ,N - 1.  Regardless 
of  the  position  of  atoms,  all  pump  pairs  with  frequency  dif¬ 
ference  nd  contribute  to  the  spin  interaction  Jn.  Considering 
first  the  XY  terms,  and  according  to  Eq.  7,  we  demand  the 
following: 


N-n-l 


jn*j  ]T  xax*+n  =— 


Jo 


a— 0 


sin2(mr/A/y 


[14] 


where  J  is  the  GM  photon  coupling  rate  (Eq.  8)  that  we  will 
assume  to  be  a  constant  for  the  simplicity  of  discussions.  This 
requires  that  the  physical  size  of  the  spin  chain  be  small 
compared  with  the  decay  length  of  J.  That  is,  Nd  <C  £  where 
d  is  the  atomic  separation.  It  is  then  straightforward  to  find 
the  required  pump  amplitudes  Qa  (or  equivalently  Xa )  by 
solving  Eq.  14  for  all  n.  Notice  that  the  system  of  equations 
Eq.  14  is  overdetermined,  and  therefore  one  can  find  several 
solutions  of  it.  However,  we  choose  the  solution  that  mini¬ 
mizes  the  total  intensity  |^«|2-  Fig-  5  shows  that  the  total  in¬ 
tensity  converges  to  a  constant  value  for  large  N,  as  a  result  of 
decreasing  sideband  amplitudes  for  decreasing  1/r2  interaction 
strengths.  This  is  confirmed  in  Fig.  5  as  we  see  the  growth  of 
the  ratio  between  maximum  and  minimum  sideband  amplitudes 
when  N  increases.  The  same  external  pump  configuration  can  also 
be  used  to  induce  the  ZZ  terms  by  applying  stroboscopic  proce¬ 
dures  as  discussed  in  the  previous  section. 

Engineering  Spin  Hamiltonians  for  2D  Systems:  Topological 
and  Frustrated  Hamiltonians 

In  the  following,  we  discuss  specific  examples  for  engineering  2D 
spin  Hamiltonians  that  are  topologically  nontrivial.  In  particular, 
we  discuss  two  chiral-flux  lattice  models  that  require  long-range 


Nt 

Hxy 

Hzz 

Hxy 

Hzz 

Hxy 

Hzz 

Time  t 

Fig.  4.  Scheme  for  generating  an  XXZ  spin  Hamiltonian  using  a  stroboscopic 
evolution.  The  scheme  contains  periodic  applications  of  a  multifrequency  Raman 
field  to  induce  the  HXy  interaction  (in  green),  two  fast  microwave  pulses  (or 
optical  two-photon  transition)  forming  Hmw  that  uniformly  rotate  the  spin  basis 
{|9)rn  |s)„}  ♦+  {(|g)„  +  |s)„)/n/2,  (Ig)„  -  |s)„)/v5}  back  and  forth  (in  blue),  and  a 
butterfly-like  pumping  scheme  that  applies  HZz  in  the  rotated  basis. 
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hopping  terms  to  engineer  single  particle  flat-bands  with  nonzero 
Chern  numbers,  which  are  key  ingredients  to  realizing  fractional 
quantum  Hall  effects  (FQHEs)  without  Landau  levels  (58,  59). 

In  recent  years,  the  field  of  ultracold  atoms  has  made  re¬ 
markable  progress  in  engineering  topological  quantum  matter. 
An  artificial  gauge  field  (67)  has  been  realized  using  cold  atoms 
loaded  into  shaken  optical  lattices  (68,  69)  as  well  as  in  lattices 
with  laser-induced  tunneling  (55,  54,  70-72).  Various  topolog¬ 
ical  models,  including  Haldane’s  honeycomb  lattice  (73,  74), 
have  been  successfully  implemented.  Berry  curvature  and  to¬ 
pological  invariants  such  as  the  Chern  number  (74-80)  can  be 
measured.  Chiral  edge  currents  in  synthetic  quantum  Hall  lat¬ 
tices  are  also  observed  (81,  82).  Most  of  the  demonstrations 
so  far  focus  on  probing  topological  band  structures  and  single¬ 
particle  physics.  Realizing  strongly  interacting  topological 
phases  such  as  FQH  states,  however,  still  remains  elusive.  This 
in  part  is  due  to  limited  topological  bandwidth-to-gap  ratio,  but 
a  number  of  improved  schemes  (e.g.,  refs.  83  and  84)  have 
been  proposed. 

Coupling  cold  atoms  to  mobile  PCW  photons  also  allows  to¬ 
pological  band  engineering  and  band  flattening.  Moreover,  the 
pseudo  spin- 1/2  system  already  interacts  like  hard-core  bosons 
because  individual  atoms  that  participate  in  the  spin-exchange 
process  cannot  be  doubly  excited.  With  the  addition  of  tunable 
long-range  ZZ  interactions,  we  can  readily  build  many-body 
systems  that  should  exhibit,  for  example,  FQH  and  supersolid 
phases  (85),  providing  a  powerful  route  toward  realizing  strongly 
interacting  topological  phases. 

Chiral-Flux  Square  Lattice  Model.  The  first  example  discussed  here 
can  be  mapped  to  a  topological  flat-band  model  similarly  de¬ 
scribed  in  refs.  58,  59,  and  85.  The  topological  spin  Hamiltonian 
is  written  as  follows: 


#fiat  =Hq+H' 

^2  °m°n+  h.C 

(m,n)  (( m,n ))  [15] 

H'  =  t3  ^2 

(((m,n))) 

in  which  we  define  0^=0^  and  am=o^s\  (.)  denotes  nearest 
neighbors  (NN),  and  t\  is  the  coupling  coefficient,  ((.))  [(((.)))] 
denotes  [next-]next-nearest  neighbors  (NNN)  [and  NNNN,  re¬ 
spectively]  with  t3  fe]  being  the  respective  coupling  coefficients. 
The  NN  coupling  phases  </>mn  =  ±<p  are  staggered  across  lattice 
sites,  where  the  phase  factor  </>  is  the  one  that  breaks  time  reversal 
symmetry  for  cf)^0,nn  (with  n  e  Z).  Spin  exchange  between  next- 
nearest  neighbors  (NNN)  has  real  coefficients  ±t2  with  alternating 
sign  along  the  lattice  checkerboard  (Fig.  6).  One  can  show  that 


N 

Fig.  5.  Sideband  amplitude  for  Haldane-Shastry  model:  total  intensity 
(black)  J2\xa\2  and  maximum/minimum  ratio  (red)  of  sideband  amplitudes 

\Xa\  as  a  function  of  N. 


Fig.  6.  Engineering  a  chiral-flux  square  lattice.  (A)  Two  sublattices  ( nx  +  ny 
odd  or  even)  are  marked  by  blue  and  red  circles,  respectively.  Solid  lines 
mark  the  NN  hopping  with  phase  gain  cp  (arbitrarily  tuned)  along  the  di¬ 
rection  of  the  arrows.  Dashed  (dotted)  lines  mark  the  NNN  hopping  terms 
(coefficients  ±£2)-  NNNN  long-range  hopping  along  curved  lines  are  included 
to  assist  band  flattening.  Filled  arrows  indicate  the  propagation  of  pump 
electric  fields  sy  and  sXl  respectively;  see  text.  ( B )  Resulting  two-band  struc¬ 
ture  with  (t2,  t3l  (p)  =  (ti  /a/2,  ti  /4a/6,  n/A). 


already  Hq  has  a  small  bandwidth  with  nontrivial  Chern  number 
that,  choosing  t2—t\j\fl  and  (f)  =  n/ 4,  results  in  a  simple  band 
dispersion  Lo(k)=  ±v^2*i  a/3  +  cos (kx  4-  ky)cos(kx  - ky).  Adding 
H'  to  H{)  with,  for  example,  t3  =  1  /4v6*i  allows  us  to  engineer  an 
even  flatter  lower  band  whose  bandwidth  is  ~  1  %  of  the  band  gap. 

We  can  use  an  array  of  atoms  trapped  within  a  2D  PCW,  as 
in  Fig.  IB,  to  engineer  the  Hamiltonian  i/fiat  of  Eq.  15.  For 
simplicity,  we  assume  that  there  is  one  atom  per  site  although 
this  is  not  a  fundamental  assumption.8  As  shown  in  Fig.  6 A,  we 
need  to  engineer  spin  exchange  in  four  different  directions, 
namely,  x,y,x±y.  We  first  introduce  linear  Zeeman  shifts  by 
properly  choosing  a  magnetic  field  gradient  VB  (SI Appendix  B: 
Proper  Choice  of  Ground- State  Energy  Shifts  in  2D  Models)  such 
that  8 a  =  \pBS7B  ■  Ara\,  where  ptB  is  the  magnetic  moment,  Ara 
are  vectors  associated  with  the  directions  of  spin  exchange: 


{Ara}„=x,>aw.  =  {dx,dy,d(x+y),d(x-y)\,  and  d  is  the  lattice 
constant.  To  activate  spin  exchange  along  these  directions  while 
suppressing  all  other  processes,  we  consider  a  simplest  case  by 
applying  a  strong  pump  field  of  amplitude  £2o  (frequency  coi)  to 
pair  with  sidebands  |£2a|  <C  |£2o|  of  detunings  cba  =  8a  to  satisfy  the 
resonant  conditions.  To  generate  the  desired  chiral-flux  lattice,  we 
need  to  carefully  consider  the  propagation  phases  ko  •  rn  (ka  •  rn) 
of  the  pump  field  (and  sidebands),  where  rn  =  d(nxx  +  nyy)  is  the 
site  coordinate  and  nxy  e  Z.  In  the  following,  we  pick  k  =  ka  =  n/d. 

We  can  generate  the  couplings  in  7/flat,  term  by  term,  as  follows. 
Staggered  NN  coupling  along  Ara=Xjy.  We  consider  the  strong  pump 
field  to  be  propagating  along y,  that  is,  Xo(rn)  =  (\Q.o\/2A)e~lny!r. 
At  the  NN  site  rm  =  rn  +  Arx,  it  can  pair  with  an  auxiliary 
sideband  of  detuning  cbx  =  8X  =  \pBV B  -  Arx\  with  Xx(rm)  = 
(|Qi  |/2 A)  [e~inyn  -  i^e-linx+1)^]  to  generate  coupling  along  Ar*.  The 
sideband  is  formed  by  two  field  components  in  sy(t)  and  sx(t),  propa¬ 
gating  along  y  and  x,  respectively  (Fig.  6),  with  an  amplitude  ratio  of 
£  and  with  an  initial  n /2  phase  difference.  These  two  fields  are  used  to 
independently  control  real  and  imaginary  parts  of  the  spin-exchange 
coefficients.  Using  Eq.  8  under  the  condition  |£2o|  >  |£2«|,  the  cou¬ 
pling  rate  along  Axx  is  as  follows: 


Jm,n=JX{ vX;=h 


Vi+c1 


=  t\e 


,±i(p 


[16] 


where  t\  =J  \Xq  \  \Xx  \  +  This  results  in  the  staggered  phase 

pattern  with  tunable  </>  =  tan-1  f .  The  NN  coupling  along  Ary  can 


§ln  principle,  exact  physical  separations  between  trapped  atoms  do  not  play  a  significant 
role  with  photon-mediated  long-range  interactions.  One  may  also  engineer  the  spin 
Hamiltonian  based  on  atoms  sparsely  trapped  along  a  photonic  crystal,  even  without 
specific  ordering.  It  is  only  necessary  to  map  the  underlying  symmetry  and  dimensionality 
of  the  desired  spin  Hamiltonian  onto  the  physical  system. 
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Fig.  7.  Engineering  a  honeycomb-equivalent  topological  brick  wall  lattice. 
(A)  Unit  cell  of  a  honeycomb  lattice.  Solid  lines  mark  the  NN  hopping. 
Dashed  lines  mark  the  NNN  hopping  with  phase  gain  $  along  the  direction 
of  the  arrows.  ( B )  Unit  cell  of  a  brick  wall  lattice.  Solid  lines  indicate  the  NN 
hopping  as  in  A.  NNNN  hopping  (curved  dashed  lines)  and  NNN  hopping 
(diagonal  dashed  lines)  correspond  to  the  complex  NNN  hopping  in  A, 
making  the  two  models  topologically  equivalent.  (C)  Brick  wall  lattice.  Filled 
arrows  illustrate  the  pump  electric  fields.  (D)  Band  structure  of  the  brick  wall 
lattice,  plotted  with  cos^  =  3^3/43  (58). 


be  introduced  via  another  sideband  with  detuning  cby  =  8y  = 
\pBVB- Ary\ and  Xy(rn  +  Ary)=  -(|£2i|/2A)  [e~liny+1)n  +i£e~irix!r]. 
NNN  coupling  along  Ara=xy,xy*.  The  sign  of  the  coefficient  depends 
on  the  sublattices.  To  engineer  these  couplings,  we  use  two 
sidebands  formed  by  field  components  in  sx(t),  with  detunings 
8a=xyjy'  and  Xa(rn  +  Arw«)  =  ±(|Q2|/2A)e_‘>'^+1)  at  NNN  sites. 
After  pairing  with  the  pump  field  Xq  at  site  rn,  the  resulting  ex¬ 
change  coefficients  are  Jm/l  =JX oX*=xy^y*  =  +t2(-l)nx~ny,  forming 
the  required  pattern  with  ^  =/|^o||^|- 

NNNN  coupling  along  2Ara=x,y.  We  use  two  sidebands  X^y  = 
|Q3|e_WI”y/2A,  propagating  along  y  with  detunings  28a=xy,  to  in¬ 
troduce  the  real  coupling  coefficient  t$  =J\Xo\\X2x\. 

Summing  up,  all  of  the  components  in  the  Raman  field  can  be 
introduced  by  merely  two  pump  beams  propagating  along  x  and 
y  directions,  respectively.  In  SI  Appendix  C:  Pump  Field  Config¬ 
urations  for  Engineering  a  Chiral-Flux  Square  Lattice  Model ,  we 
explicitly  write  down  the  time-dependent  electric  field  that 
contains  all  of  the  sidebands. 

We  note  that  it  is  also  possible  to  simultaneously  introduce  both 
blue-detuned  (<5a>0)  and  red-detuned  ( 8_a  =  -8a )  sidebands  in 
the  Raman  field  to  control  the  same  spin-exchange  term.  That  is, 
Jm,n  =J[Xo(rn)X*(rm)  +Xq  (rm)X_a(r„)],  which  has  contributions 
from  Xa  and  X-a  of  blue  and  red  sidebands,  respectively.  Arranging 
both  sidebands  with  equal  amplitudes  lead  to  equal  contribu¬ 
tions  in  the  engineered  coupling  coefficient.  This  corresponds 
to  applying  amplitude  modulations  in  the  pump  electric  field. 
In  real  experiments,  amplitude  modulation  can  be  achieved 
by,  for  example,  the  combination  of  acoustic-optical  modulators, 
and  optical  IQ-modulators. 


NNNN  (along  2Ary)  couplings  with  the  same  strength  and 
with  a  coupling  phase  (f)mn  =  ±<p,  which  alternates  sign  across 
two  sublattices.  Thus,  our  target  Hamiltonian  is  given  by  the 
following: 


H  —  h  {pln^n  +h.C.)  +t2  (e1^™ <7^(7 n  +h.c.),  [17] 

(m,n)  {m,n} 


where  (. )  denotes  NN  pairs  in  the  brick  wall  configuration  (Fig. 
7)  and  t\  is  the  coupling  coefficient.  Note  that,  for  simplicity,  we 
discuss  a  special  case  where  all  NN-coupling  coefficients  from  a 
brick  wall  vertex  are  identical.  The  second  summation  in  Eq.  17 
runs  over  both  NNN  and  NNNN  pairs  with  identical  coupling 
coefficient  ^  and  alternating  phase  (j)mn  =  ±(p  (Fig.  7). 

As  in  the  previous  case,  we  use  a  strong  pump  field  (propa¬ 
gating  along  y),  as  well  as  several  other  weak  sidebands  to  gen¬ 
erate  all  necessary  spin-exchange  terms.  Detailed  descriptions  on 
engineering  individual  terms  can  be  found  in  SI  Appendix  D: 
Pump  Field  Configurations  for  Engineering  a  Topological  Spin 
Model  in  a  Brick  Wall  Lattice.  The  most  important  ingredient, 
discussed  here,  is  that  we  can  generate  checkerboard-like  NN 
coupling  (along x),  with /m>„  =JXqX*  =  {t\  /2)  [1  -  (~l)nx~ny].  This  is 
achieved  by  using  a  sideband  of  detuning  8X  and  amplitude 
Xx  =  (| £2 1/4 A)  [e~inyn  +  ^e-i{nx+i)n]  at  position  rm  =  rn  +  Arx,  formed 
by  two  fields  propagating  along  y  and  x,  respectively.  If  both 
fields  have  the  same  amplitude  (f  =  1),  they  either  add  up  or 
cancel  completely  depending  on  whether  nx  -  ny  is  odd  or  even. 
If  one  applies  the  same  trick  toward  NN  coupling  along  y,  but 
with  f^l,  the  coupling  amplitude  modulates  spatially  in  a 
checkerboard  pattern.  Essentially,  all  three  NN  terms  around  a 
brick  wall  vertex  can  be  independently  controlled,  opening  up 
further  possibilities  to  engineer,  for  example,  Kitaev’s  honey¬ 
comb  lattice  model  (87,  88). 

For  physical  implementations,  again  only  two  pump  beams  can 
introduce  all  components  required  in  the  Raman  field,  which  is 
very  similar  to  the  previous  case.  We  stress  that,  by  merely 
changing  the  way  the  Raman  field  is  modulated,  one  can  dy¬ 
namically  adjust  the  engineered  spin  Hamiltonians  and  even 
the  topology,  as  we  compare  both  cases.  This  is  a  unique  feature 
enabled  by  our  capability  to  fully  engineer  long-range  spin 
interactions. 

Moreover,  many  of  the  tricks  discussed  above  can  also  be 
implemented  in  ID  PCWs.  It  is  even  possible  to  engineer  a 
topological  ID  spin  chain,  by  exploiting  long-range  interac¬ 
tions  to  map  out  nontrivial  connection  between  spins.  For 
example,  our  method  can  readily  serve  as  an  realistic  ap¬ 
proach  to  realize  a  topological  ID  spin  chain  as  recently 
proposed  in  ref.  89. 

XXZ  Spin  Hamiltonian  with  Tunable  Interaction  1  /r*1.  In  the  last  ex¬ 
ample,  we  highlight  the  possibility  of  engineering  a  large  class  of 
XXZ  spin  Hamiltonians,  which  were  studied  extensively  in  the 
literature  because  of  the  emergence  of  frustration  related  phe¬ 
nomena  (60,  90-96)  and  their  intriguing  nonequilibrium  dynamics 
(45-50).  An  XXZ  Hamiltonian  is  typically  written  as  follows: 

Hxxz  =  ~B  <  +  E  [«>s(0K<4  +  sin(<?)  (dndm  +  «,)] , 

n  n<m '  n’m 

[18] 


"Honeycomb"-Equivalent  Topological  Lattice  Model.  To  further 
demonstrate  the  flexibility  of  the  proposed  platform,  we  create 
Haldane’s  honeycomb  model  (73)  via  a  topologically  equivalent 
brick  wall  lattice  (74,  86).  Here,  we  engineer  the  brick  wall 
configuration  using  the  identical  atom-PCW  platform  discussed 
in  the  previous  example.  Mapping  between  the  two  models  is 
illustrated  in  Fig.  7  A  and  B ,  which  contains  the  following  two 
nontrivial  steps:  (/)  generating  a  checkerboard-like  NN-exchange 
pattern  in  the  x  direction;  (//)  obtaining  NNN  (along  Ar^,*)  and 


where  an  effective  magnetic  field  B  controls  the  number  of  excita¬ 
tions,  rm,m  =  \rn  -  rm|,  and  the  parameter  6  determines  the  relative 
strength  between  the  ZZ  and  XY  interactions.  This  class  of  spin 
models  has  been  previously  studied,  but  mostly  restricted  to  nearest 
neighbors  (90-93)  or  dipolar  ( rj  =  3)  interactions  (60,  94,  95). 

In  our  setup,  one  can  simulate  XXZ  models  with  arbitrary  rj  by 
first  introducing  unique  ground-state  energy  shifts  at  each  of  the 
separation  rn  -  rm,  and  then  applying  a  strong  pump  field  of  am¬ 
plitude  £2q  together  with  Nd  auxiliary  fields  Qa  of  different  detunings 
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Fig.  8.  Mean  magnetization  M/N  for  a  system  with  N=  16atoms  in  a  square 
lattice,  restricted  to  A/exc<8  excitations  and  rj=  1  (A),  2  ( B ),  3  (C),  and  NN 
couplings  (D). 


of  the  many-body  states  may  be  probed  via  guided  photons  in  an¬ 
other  propagation  mode  along  the  PCW  (98). 

Limitations  and  Error  Analysis.  Until  now,  we  have  mainly  focused 
on  how  to  engineer  Hq  in  an  ideal  situation.  We  neglected  spon¬ 
taneous  emission  or  GM  photon  losses  and  considered  that  the 
energy  gradient  (or  8,  the  ground-state  energy  difference  between 
nearest  neighboring  atoms)  can  be  made  very  large  compared 
with  the  interaction  energy  scales  that  we  want  to  simulate 
(\8\  >  \Jm,m+i  !)•  Because  the  effect  of  finite  cooperativities  was 
considered  in  detail  in  refs.  36  and  37,  and  their  conclusions 
translate  immediately  to  our  extension  to  multifrequency  pumps, 
in  this  work  we  mainly  focus  on  the  effect  of  finite  8.  In  addition, 
we  also  discuss  the  effects  of  AC  Stark  shifts  as  in  Eq.  9,  and  its 
error  contributions,  together  with  other  possible  error  sources. 

Corrections  Introduced  from  Higher  Harmonics:  A  Floquet  Analysis. 

We  discuss  errors  and  the  associated  error  reduction  scheme 
following  a  Floquet  analysis  with  multifrequency  driving  (99, 
100),  applicable  mainly  to  ID  models.  Including  all  of  the  time- 
dependent  terms  in  a  multifrequency  pumping  scheme,  we  have 
(Eq.  5)  H(t)  =  Y  HpeipSt,  where  Hp  represents  the  part  that  oscil¬ 
lates  at  frequency  pS.  This  Hamiltonian  has  a  period  T  =  2ji/8.  It 
can  be  shown  that  at  integer  multiples  of  T,  the  observed  system 
should  behave  as  if  it  is  evolving  under  an  effective  Hamiltonian#: 


to  introduce  spin  interactions  at  each  separations.11  Moreover,  the 
parameter  6  that  determines  the  ratio  between  ZZ  and  XY  in¬ 
teraction  can  be  controlled  by  using  different  pump  intensities  in  the 
stroboscopic  steps  (SI  Appendix  E:  PCW  and  Pump  Field  Configu¬ 
rations  for  Engineering  an  XXZ  Spin  Hamiltonian  with  1/ rn 
Interaction). 

To  illustrate  physics  that  can  emerge  in  the  first  experimental 
setups  with  only  a  few  atoms,  we  study  the  total  magnetization  of  a 
small  square  lattice  of  ns  x  ns  (=N)  16  atomic  spins.  We  apply  exact 
diagonalization  restricting  to  Nexc  <  8  excitations  for  N  =  16  spins 
and  cover  one-half  of  the  phase  diagram  with  B  >  0.  In  Fig.  8,  we 
explore  the  mean  magnetization  of  the  system  M/N  =  \  Y?  (df)  /N  as 
a  function  of  B  and  6  for  rj  =  1  (A),  2  ( B ),  3  (C),  and  NN  cou¬ 
plings  (D).  At  6  =  0,  the  system  behaves  classically  showing  the  so- 
called  “devil  staircase”  (97)  of  insulating  states  with  different 
rational  filling  factors  and  crystalline  structures.  As  already  explored 
in  ref.  95  for  ID  dipolar  systems  (rj  =  3),  the  presence  of  long-range 
interactions,  compared  with  nearest-neighbor  models,  lead  to 
stronger  frustration  effects.  This  manifests  in  the  magnetization 
diagram  with  an  asymmetry  between  6  $  0.  In  Fig.  8,  we  show  that 
longer-range  interactions  lead  to  even  higher  degree  of  asymmetry. 
Moreover,  in  refs.  60  and  95,  it  was  discussed  that  long-range 
coupling  leads  to  the  formation  of  supersolid  phases,  in  which 
crystalline  structure  and  long-range  order  coexist.  These  may  be 
even  more  favored  by  longer-range  interactions.  Full  character¬ 
ization  of  the  phases  diagram  is,  however,  beyond  the  scope  of  this 
paper  and  will  be  discussed  elsewhere. 

Another  especially  interesting  arena  is  the  behavior  of  strongly 
long-range  interacting  systems  ( rj  smaller  than  the  lattice  di¬ 
mension  D )  under  nonequilibrium  dynamics.  It  has  recently  been 
predicted  to  yield  “instantaneous”  transmission  of  correlations 
after  a  local  quench  (45,  47,  48,  96),  breaking  the  so-called  Fieb- 
Robinson  bound. 

Finally,  it  is  interesting  to  point  out  that  magnetization  can  be 
measured  by  first  freezing  the  interaction  (via  shutting  off  the  pump 
lasers)  followed  by  atom  number  counting  using  state-dependent 
fluorescence  imaging.  Coherence  and  off-diagonal  long-range  orders 


#eff,l 


[HP,H-p] 

P 


1  v[[j Hp,Ho],H-p\  +  [[H-p,Ho\,Hp\ 
2  S2^  P2 


[19] 


This  means  that  the  leading  error  in  our  simple  scheme  would  be 
on  the  order  of  J2/8,  where  J  is  the  simulated  interaction 
strength.  However,  we  note  that  if  Hp  =  ±H-P,  the  leading  error 

term  Y[Hp>H-p\/ (P&)  should  vanish.  In  other  words,  first-order 

error  vanishes  if  Hp  is  either  symmetric  or  antisymmetric  under  a 
time  reversal  operation  ST.  Although  the  original  Hamiltonian 
H(t)  does  not  necessarily  possess  such  symmetry,  it  is  possible 
to  introduce  a  two-step  periodic  operation  instep  =  {H,  SEH, 
. . . }  to  cancel  the  first-order  error  while  keeping  the 
time-independent  part  //2step,o=#o  identical.  This  results  in  an 
effective  Hamiltonian  in  the  Floquet  picture: 

Hef f,2  =H0  +  Hqxx2  «H0+I  Y.  Hf  UElEAA,  [20] 

8  z '  D 

p  1 

where  Hp  is  the  (operator)  Fourier  coefficient  of  the  two-step 
Hamiltonian  and  the  leading  error  reduces  to  the  order  of  J3  /82. 

To  achieve  the  time  reversal  operation,  we  must  reverse  the 
phase  of  the  driving  lasers,  as  well  as  the  sign  of  the  energy 
offsets  between  the  atoms.  Specifically,  we  can  engineer  a  peri¬ 
odic  two-step  Hamiltonian  by  first  making  the  system  evolve 
under  presumed  Ho  (along  with  other  time-dependent  terms)  for 
a  time  interval  T,  and  then,  for  the  next  time  interval  T,  we  flip 
the  sign  of  the  energy  gradient,  followed  by  reversing  the  prop¬ 
agation  direction  of  the  Raman  fields  such  thatXa  -+X*  in  Eq.  5. 
As  a  result,  all  of  the  time-dependent  Hamiltonians  Hp ,  \/p  ±  0, 
become  H_p  in  the  second  step,  resulting  in  Hp  =  (-1  ffH~p  re¬ 
quired  for  error  reduction;  whereas  the  time-independent 
Hamiltonian  Ho  remains  identical  in  the  two-step  Hamiltonian. 
See  SI  Appendix  F:  Error  Reduction  and  Analysis  for  more  discussions. 


^To  simulate  a  square  lattice  of  nsxns  ( =N)  atomic  spins,  we  find  that  the  number  of 
different  distances  grows  as  A/d  =  (ns(ns  +  1)-2)/2f  which  is  linearly  proportional  to  the 
number  of  atoms  Ndoc.N. 


#When  the  measurement  time  is  incommensurate  with  period  T,  small-amplitude  and 
fast-oscillating  spin-dynamics  due  to  time-dependent  terms  in  Eq.  5  manifest  as  extra 
errors;  see  the  discussion  about  micromotion  in  refs.  99  and  100. 
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Fig.  9.  (/A)  Comparison  of  ground-state  energy  error  |(£0 -£,-)/£o|  and  ( B ) 

ground-state  overlap  |('P0|'P/)|  as  a  function  of  N  for  Hamiltonians  Heff(1 
(black)  and  He ffi2  (red)  with  detuning  6jJ\  =40. 


Numerical  Analysis  on  the  Haldane-Shastry  Spin  Chain.  We  now 

analyze  numerically  and  discuss  error  on  one  particular  example. 
For  numerical  simplicity,  we  choose  the  Haldane-Shastry  model 
as  its  ID  character  makes  it  numerically  more  accessible.  How¬ 
ever,  the  conclusions  regarding  the  estimation  of  errors  can  be 
mostly  extended  to  other  models.  As  we  have  shown  in  Eq.  13, 
the  Haldane-Shastry  Hamiltonian  is  composed  by  an  XY  term 
plus  a  ZZ  term  that  we  can  simulate  stroboscopically.  As  we 
already  analyzed  the  Trotter  error  due  to  the  stroboscopic  evo¬ 
lution,  here  we  focus  on  the  XY  part  of  the  Hamiltonian,  which 
reads  as  follows: 


Jo 


sin  (nn/N) 


[21] 


where  Jo  =Jn2/N2.  Following  the  prescribed  engineering  steps, 
the  total  time-dependent  Hamiltonian  resulting  from  multiple 
sidebands  can  be  written  as  H(t)  =  Y/  Hpeip5t,  with  the  following: 

p 

N  N-m 

hp=EE  (4(P)«+n  +4A)<«)’  [22] 

m— 1  n— 1 

and  we  have  defined  J„,{p)  =  Y,aJ)LXxp8n-p,/i-a-  Here,  Xa  are 
fixed  such  that  Ho  =Hus,xy 

To  illustrate  the  effect  of  error  cancellations,  we  consider  first 
a  scenario  where  we  directly  apply  Eq.  22.  To  leading  order,  the 
effective  Hamiltonian  is  //eff,i  as  in  Eq.  19.  We  then  analyze  the 
two-step  driving,  using  the  effective  Hamiltonian  //eff,2  in  Eq.  20, 
with  Hp  given  by  Eqs.  S45  and  S46. 

We  calculate  the  ground-state  energies  and  eigenvectors  of 
Ho ,  He ff;i,  and  He ff;2,  which  we  denote  as  E0)  1)2  and  |¥0jij2),  re¬ 
spectively,  for  different  number  of  atoms  and  different  ratios  of 
S/Jo.  The  results  are  shown  in  Figs.  9  and  10.  In  panels  A ,  we 
show  the  error  in  absolute  value  with  respect  to  the  ideal 
Hamiltonian  Ho.  Interestingly,  due  to  particular  structure  of  |To) 
and  Hp ,  one  can  show  that  (T0|[F/p,//_;7]|To)  and  the  first- 
order  correction  to  the  energy  vanishes.  This  is  confirmed  in  Fig. 
10 L4,  where  we  found  that  the  error  actually  scales  with  1  /S2. 

Moreover,  it  is  also  enlightening  to  compare  the  overlap  of  the 
ground  states  as  shown  in  Figs.  9 B  and  10 B.  We  only  compute 


the  even-atom  number  configuration  as  the  odd  ones  are  de¬ 
generate  and  therefore  the  ground  state  is  not  uniquely  defined. 
We  see  that  the  ground-state  overlap  of  //eff,2  is  several  orders 
of  magnitude  better  than  the  one  with  Moreover,  its 

dependence  on  S  is  better  than  the  l/<52  expectation. 

The  Role  of  Time-Dependent  Stark  Shifts  in  the  Error  Analysis.  In  the 

previous  discussions,  we  have  dropped  the  contribution  of  the 
time-dependent  Stark  shifts: 


mp-l 

Hac(t)=-j2J2* 

n  a>j3 


2A 


piVajt 


[23] 


where  6)ajj  =  cba~  d>p-  In  SI  Appendix  F:  Error  Reduction  and 
Analysis ,  we  discuss  its  role  in  the  effective  Hamiltonian,  using 
the  Floquet  error  analysis.  To  summarize,  we  evaluated  the  error 
in  the  two-step  driving  scheme  in  various  configurations. 

Generic  Hamiltonians  with  translational  invariance.  By  translational 
invariance,  we  mean  that  there  are  no  site-dependent  spin  in¬ 
teractions,  and  the  spin-exchange  coefficients  remain  identical  as 
we  offset  the  spin  index  by  one  or  more.  This  means  that  all 
components  in  the  pump  field  should  drive  the  system  with 
uniform  optical  phases  as  in  the  Haldane-Shastry  model  dis¬ 
cussed  above.  The  error  by  Hac(t)  averages  out  to  zero  in  the 
Floquet  picture.  In  the  butterfly  scheme,  however,  both  | g)  and 
| s)  states  are  pumped  and  they  may  be  shifted  differently.  This 
leads  to  slight  modifications  in  the  engineered  XX  and  YY  terms. 
Models  containing  sublattices.  For  topological  models  that  contain 
sublattices,  as  in  our  examples,  the  pump  fields  are  not  perfectly 
transverse  and  Stark  shifts  are  site  dependent,  resulting  in  non¬ 
vanishing  error.  For  realistic  PCW  realizations,  one  should  set 
moderate  pump  detuning  J / A  >  0(1)  such  that  leading  error 
contribution  will  be  <J3/S2,  and  for  <5»/,  the  Stark  shift 
terms  may  be  ignored. 

Stark  shift-dominated  regime.  It  may  be  possible  that  our  sublattice 
models  be  purposely  driven  with  large-amplitude  pumps  such 
that  |£2|2/A  >  S.  Stark  shift  contributions  would  become  important 
in  the  resulting  spin  dynamics.  However,  if  we  choose  a  large  pump 


S/Ji 


Fig.  10.  (A)  Comparison  of  ground-state  energy  error  |(£0-£,)/£o|  and  (£) 
ground-state  overlap  |('P0|'P/)|  as  a  function  of  for  Hamiltonians  Heffi1 
(black)  and  Hef f  2  (red)  for  A/  =  12  atoms. 
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detuning  A  >7,  the  dominant  error  contribution  can  in  fact  be 
written  in  the  following  simple  form: 

He„a  *  ^Am,n  (jm/logo?g  +h.c .) ,  [24] 

m,n 

where  Am,n  is  a  site-dependent  amplitude.  In  a  special  case  that 
only  two  sublattices  are  present,  as  in  our  examples,  we  note  that 
Am^n  may  only  depend  on  the  distance  and  is  site  independent. 
This  “error”  term  would  then  uniformly  modify  the  XY  coupling 
strengths  to  a  new  value: 

J m,n  =  (l  m,n)Jm,ri'  [25] 

The  next  leading  order  errors  are  a  factor  of  ~7/A  smaller  than 
this  leading  Stark  shift  contribution,  suggesting  we  can  always 
increase  the  detuning  A,  while  keeping  |£2|  /A  constant,  to  reduce 
the  error  contribution. 

Other  Error  Sources  and  Heating  Effects.  Apart  from  errors  arising 
from  multifrequency  driving,  there  are  other  common  error 
sources  in  cold  atoms  that  we  have  not  considered  so  far,  such  as 
motional  heating.  In  the  PCW  platform,  atoms  are  tightly  con¬ 
fined  with  a  trap  depth  more  than  three  orders  of  magnitude 
larger  than  the  recoil  energy,  rendering  well-separated  motional 
bands  such  that  effects  like  interband  heating  (101)  can  be 
suppressed.  Spin-exchange  rates  in  the  PCW  platform,  however, 
can  be  adjusted  to  1  MHz  >  \Jxy’z\  >  1  kHz  so  that  the  many- 
body  time  scales  1  ms)  can  be  much  faster  than  those  asso¬ 
ciated  with  motional  heating. 

In  fact,  spin  temperature  can  be  decoupled  from  real  atomic 
temperature  while  simulating  the  spin  models.  For  example,  one 
can  polarize  atomic  spins  initially  in  a  strong  magnetic  field 
( B  >  to  approximate  a  zero-temperature  paramagnetic 

phase  (17).  The  magnetic  field  can  then  be  ramped  down  adia- 
batically  to  the  final  value  of  the  desired  spin  model.  Limitations 
to  adiabaticity  and,  therefore,  to  the  accessible  spin  temperature 
will  ultimately  be  limited  by  the  fidelity  of  the  spin-exchange  (36, 
37)  or  by  motional  heating  that  leads  to  dephasing,  whichever 
gives  a  more  stringent  bound. 

Conclusions  and  Outlook 

In  this  paper,  we  have  shown  that  atom-nanophotonic  systems 
present  appealing  platforms  to  engineer  many-body  quantum 
matter  by  using  low-dimensional  photons  to  mediate  interaction 
between  distant  atom  pairs.  We  have  shown  that,  by  introducing 
energy  gradients  in  ID  and  2D,  and  by  applying  multifrequency 


Raman  addressing  beams,  it  is  possible  to  engineer  a  large  class 
of  many-body  Hamiltonians.  In  particular,  by  carefully  arranging 
the  propagation  phases  of  Raman  beams,  it  is  possible  to  introduce 
geometric  phases  into  the  spin  system,  thereby  realizing  nontrivial 
topological  models  with  long-range  spin-spin  interactions. 

Another  appealing  feature  of  our  platform  is  the  possibility  of 
engineering  periodic  boundary  conditions,  as  explicitly  shown  in 
the  ID  Haldane-Shastry  model,  or  other  global  lattice  topology  by 
introducing  long-range  interactions  between  spins  located  at  the 
boundaries  of  a  finite  system.  Using  2D  PCWs,  for  example,  it  is 
possible  to  create  previously  unavailable  spin-lattice  geometries  such 
as  Mobius  strip,  torus,  or  lattice  models  with  singular  curvatures  such 
as  conic  geometries  (102)  that  may  lead  to  localized  topological 
states  with  potential  applications  in  quantum  computations. 

We  emphasize  that  all  of  the  pairwise- tunable  interactions  can 
be  dynamically  tuned  via,  for  example,  electro-optical  modula¬ 
tors  at  timescales  much  faster  than  that  of  characteristic  spin 
interactions.  Therefore,  the  spin  interactions  can  either  be  adi- 
abatically  adjusted  to  transform  between  spin  models  or  even  be 
suddenly  quenched  down  to  zero  by  removing  all  or  part  of  the 
Raman  coupling  beams.  We  may  monitor  spin  dynamics  with 
great  detail:  after  we  initially  prepare  the  atomic  spins  in  a 
known  state  by,  say,  individual  or  collective  microwave  address¬ 
ing,  we  can  set  the  system  to  evolve  under  a  designated  spin 
Hamiltonian,  followed  by  removing  all  of  the  interactions  to 
“freeze”  the  dynamics  for  atomic  state  detection.  Potentially,  this 
allows  for  detailed  studies  on  quantum  dynamics  of  long-range, 
strongly  interacting  spin  systems  that  are  driven  out-of-equilibrium. 
The  dynamics  may  be  even  richer  because  the  spins  are  weakly 
coupled  to  a  structured  environment  via  photon  dissipations.  We 
expect  such  a  platform  may  bring  novel  opportunities  to  the  study  of 
quantum  thermalization  in  long-range  many-body  systems,  or  for 
further  understanding  of  information  propagation  in  a  long-range 
quantum  network. 
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SI  Summary 

In  SI  Appendix  A:  Complete  Derivation  of  Final  Time-Dependent 
Hamiltonian ,  we  give  a  detailed  derivation  of  the  time-  dependent 
effective  Hamiltonian  (Eq.  4  of  the  main  text): 

N  rap- 1 

Hxr(t)=  E  [SI] 

m,  nj^m  a,  /i=0 

In  S7  Appendix  B:  Proper  Choice  of  Ground-State  Energy  Shifts  in 
2D  Models ,  we  discuss  how  to  properly  introduce  ground-state 
energy  shifts  to  engineer  generic  spin  models  in  2D.  In  SI  Appendix 
C:  Pump  Field  Configurations  for  Engineering  a  Chiral-Flux  Square 
Lattice  Model ,  SI  Appendix  D:  Pump  Field  Configurations  for  Engi¬ 
neering  a  Topological  Spin  Model  in  a  Brick  Wall  Lattice ,  and  SI 
Appendix  E:  PCW  and  Pump  Field  Configurations  for  Engineering 
an  XXZ  Spin  Hamiltonian  with  1/C  Interaction ,  we  describe  in  detail 
the  PCW  and  the  pump  field  configurations  to  engineer  a  topolog¬ 
ical  chiral-flux  lattice  model,  a  brick  wall  lattice  model,  and  an  XXZ 
model  with  1/P1  dependence,  respectively.  In  SI  Appendix  F:  Error 
Reduction  and  Analysis ,  we  discuss  in  detail  regarding  error  reduc¬ 
tion  and  analysis. 

SI  Appendix  A:  Complete  Derivation  of  Final  Time- 
Dependent  Hamiltonian 

The  PCWs  support  localized  one  or  2D  photonic  guided  modes 
(GMs),  which  can  be  described  by  a  Hamiltonian  (using  h  =  1): 

HgM  =  [S2] 

k 

where  ajk  is  the  dispersion  relation  of  the  GMs.  Neglecting  counter 
rotating  terms,  the  light-matter  Hamiltonian  can  be  written  as  follows: 

Him  —  ^^§k(y'n)^'k0^g  +  h.C.,  [S3] 

k,n 

where  gk(r«)  =gkFkTn  is  the  single-photon  coupling  constant.  The 
atomic  Hamiltonian  is  given  by  the  following: 

=  E  +  Mgf'Jgg)  ’  [S4] 


where  it  is  important  to  highlight  that  we  introduce  a  site-depen¬ 
dent  energy  in  the  hyperfine  level  | g)n  that  can  be  achieved,  for 
example,  by  introducing  a  magnetic  field  gradient  (or  a  Stark 
shift  gradient)  in  either  ID  or  2D  as  depicted  in  Fig.  1  A  and 
B  in  the  main  text.  This  site-dependent  energy,  together  with  a 
multifrequency  driving  for  \s)n  ++  \e)n  are  the  key  ingredients  of 
our  proposal.  Multifrequency  driving  with  mp  different  compo¬ 
nents  can  be  described  through  a  Hamiltonian: 

Hi  (0  =  E  <eeimL‘  +  h-cj) ,  [S5] 

where  we  have  used  the  notation  (Aab  =  \a)n  (b |,  and  col  is  the  main 
driving  frequency.  All  components  of  the  driving  field  are  em¬ 
bedded  in  Q(t),  which  can  be  written  as  follows: 

rap— 1 

0(0  s  Q.ael0<"‘,  [S6] 

a— 0 

where  cba  are  mp  different  frequency  detunings  (oba= o  =  0)  and  Qa 
the  Rabi  frequency  that  will  be  used  to  achieve  full  control  of  the 


atom-atom  interactions.  The  dynamics  of  the  system  is  described 
by  the  sum  of  all  of  the  above  Hamiltonians:  H=Ha+HoM  + 
Hdif)  +H\m. 

We  are  interested  in  the  conditions  where  |  A|  =  \coe  -  coL\  »  Q, 
such  that  the  excited  states  are  only  virtually  populated. 
To  adiabatically  eliminate  states  \e)n,  it  is  convenient  to 
work  in  a  rotating  frame  defined  by  the  transformation 
U  =  Qxp(i(J2n00^0eet  +  i'Ek®^kflk)4  which  transforms  the 
Hamiltonian  by  H  ->  UHLr  -iUdtW.  Writing  each  of  the  trans¬ 
formed  Hamiltonians,  we  have  the  following: 

Hm  -Wlm  =  +h.c., 

k  ,n 

^-^d  =  E(£y^<,+h-c)’  [ST] 

tfa-tfa  =  E(A^+cy^)’ 


while  HGm  transforms  to  zero.  Notice  that,  due  to  the  multifre¬ 
quency  driving,  it  is  not  possible  to  find  a  reference  frame  where 
the  Hamiltonian  is  time  independent.  Despite  the  time  depen¬ 
dence,  it  is  still  possible  to  adiabatically  eliminate  the  excited 
states.  For  this  purpose,  we  define  a  projector  operator  for  the 
atomic  subspace,  P  =  +  Pr°jects  out  the  excited 

states,  and  its  orthogonal  counterpart,  Q  =  ^2no^e.  Using  these 
operators,  one  can  formally  project  out  slow  and  fast  subspaces 
in  the  Schrodinger  equation: 


i^M  =  -pHF\'¥)  +PffQ|T), 
id^p-  =  QHF\'f>}+QHQ\'f>). 


[S8] 


By  using  the  fact  that  QH Q  is  actually  time  independent  and 
assuming  that  initially  there  are  no  contributions  from  the  ex¬ 
cited  states,  that  is,  Q|T(0))  =  0,  one  can  formally  integrate  Q|T) 
(by  parts),  input  the  result  into  the  equation  of  P|T),  and  obtain 
an  effective  Hamiltonian  for  the  slow  subspace: 


i^LJm  p- 

dt  \ 


QHF  P|T)  =He ffP|T).  [S9] 


The  resulting  effective  Hamiltonian  then  reads  as  the  following: 

^eff  =^eff,a  T^eff,lm  +H 


—  X!  (  0)S^n(J88^n 


Q.(t)Q.*(t) 
4A  ' 


-  + h.c. , 


[S10] 


-E 


|gk(rB)r 


k,n 


where  we  have  absorbed  some  irrelevant  phases.  The  contribution 
of  ^  =  -Ek,n(l^k(r«)|2/A)«IakO^  will  be  negligible  because  it  is 
proportional  to  the  number  of  GM  photons,  which  is  close  to  0  in 
our  situation  (36,  37).  Rewriting  the  effective  Hamiltonian  in  the 
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interaction  picture  with  respect  to  He ffja,  we  arrive  at  the  follow¬ 
ing  light-matter  Hamiltonian: 

HeCCM(t)  =  -  +  h.c.  [Sll] 

k  ,n 


coinciding  with  the  expressions  obtained  in  refs.  36  and  37  when 
\cbp  —  cogm  |  <C  | cql  —  coc  | .  Then,  depending  on  the  dimensionality 
of  the  reservoir,  we  have  the  following  (36,  37): 

/m,„,idoce-lr'-|A,  [S19] 


Note  that,  for  simplicity  in  the  derivation,  we  have  neglected  the 
contribution  of  the  time-dependent  Stark  shift  in  the  \s)n  states, 
which  is  given  by  the  following: 


2A 

[S12] 

Here,  9l[.]  indicates  real  part.  The  time-independent  contribu¬ 
tion  can  be  absorbed  into  the  energy  of  cos  without  significant 
contribution  to  the  dynamics,  whereas  the  time-dependent  terms 
will  be  averaged  out  in  the  atomic  timescales  that  we  are  inter¬ 
ested  in.  We  consider  its  possible  detrimental  effects  in  SI  Ap¬ 
pendix  F:  Error  Reduction  and  Analysis ,  where  we  analyze  the 
limitations  and  other  error  sources. 

The  relaxation  timescales  of  the  GMs  in  the  PCWs  are  typi¬ 
cality  much  faster  than  the  atomic  ones,  such  that  we  can  trace  out 
the  photonic  information  to  obtain  an  effective  master  equation 
that  describes  the  dynamics  of  the  atomic  system  through  its 
density  matrix  evolution  (103): 


8cos{t)  = 


Q(t)Q*  (t) 

4A 


mp-l 


mp-l 


=-E^-£* 

a— 0  a>P 


Jm,n,2d  [S20] 

where  =rm-rn  andXp(x)  is  the  modified  Bessel  function  of 
the  second  kind;  £  =  yjA/ 1  Axy  |  controls  the  effective  range  of  the 
interaction;  A^  =  \coc-coL  \  is  the  effective  detuning  with  respect 
to  the  band  edge,  and  A  is  the  curvature  of  the  band  (Fig.  1C). 
Notice  that  there  is  another  underlying  assumption  in  the  deri¬ 
vation,  namely,  that  the  coupling  strength  |gk|2  of  the  driven 
GMs  must  be  approximately  constant  for  all  of  the  sideband 
frequencies  col  +  oba.  If  not,  the  variation  of  |gk|2  can  be  compen¬ 
sated  by  adjusting  sideband  amplitudes  Qa.  Furthermore,  be¬ 
cause  we  focus  on  full  control  introduced  by  multifrequency 
drivings,  we  will  not  specify  the  form  of  Jmqi  and  simply  assume 
a  constant  Jm,n  =J  throughout  the  interaction  range  considered. 
However,  one  should  be  aware  that  the  length  scale  £  will  pose 
the  ultimate  limitation  of  the  range  of  the  interactions  that  we 
can  simulate. 

Because  we  have  n  (t)  =  —Tn/n  (t),  the  evolution  of  the 
density  matrix  in  Eq.  S13  is  governed  by  an  effective  XY 
Hamiltonian: 


4A2 

m,n^m  a,  ft 


°  m,n'- 


“6 


[S13] 


where  the  time-dependent  coefficients  are  given  by  the  following: 


-j£ 


F m,n{f) —  /  dsf\mne 


,i (<%m -lOgj, )<£-> (°>h +<%»! -Oh (*  S) 


4A2 


[S14] 


with  /k,m,«  =  Ek^k|2e'k(r”  Expanding  Q.{t)ST{t-s)  = 

we  find  that 


MM 


f  (A  —  _ - _ H-P  „  ei{og/n-as^+mlt-af)t 

1  m,n  \*J  —  /  J  4  A2  1  ’ 


a,p 


[S15] 


where  Tm/lfiO0  is  the  time-independent  contribution  that  can  be 
written  as  follows: 


Interestingly,  if  we  choose  a)q,n  -  a>q,m  ±  0,  we  can  control  the 
resonant  processes  by  adjusting  the  laser  frequencies  cba.  In  par¬ 
ticular,  two  atoms  n  and  m,  will  be  interacting  through  a  resonant 
process  with  rate  £la£l*pJm,n/( 4 A2)  when  the  resonant  condition, 

—  & g,n  —  dpj  —  (ba,  [S22] 

is  satisfied.  The  intuitive  picture  is  depicted  in  Fig.  7A  in  the  main 
text:  the  atom  n  scatters  from  sideband  a  a  photon  with  energy 
coL  +  d?a-  cog,n  into  GMs.  When  this  photon  propagates  to  atom 
m ,  it  will  only  be  absorbed  via  a  sideband  p  that  satisfies 
l  +  &  a  ~  <£g,n  =  eoL  +  cGp-  cog,m ,  with  the  rest  of  the  sidebands 
being  off-resonant;  Fig.  2 B  depicts  the  reversed  process.  There¬ 
fore,  the  Hamiltonian  Hxy(t)  can  be  separated  into  time-inde¬ 
pendent  (on-resonant)  and  time-independent  (off-resonant) 
contributions:  H^[t)  =Hxy$  +HX y(t),  where  7/Xy,o  is  an  XY  spin 
Hamiltonian, 


L  m,n,/3, oo 


-f 


dsf k^ 


e  -i  {o)k+o)g,m  -co L -obp)s 


[S16] 


N 

77xy,0  —  ^  Jm,nd£sdsg’  [S23] 

m,n^m 


Using  that  J0°°  eaTdz  =  n8(x)  -i?(l/x),  and  assuming  that  we  are 
working  in  the  regime  within  the  band  gap  (Fig.  ID;  assuming  an 
infinite  structure),  that  is,  Am^  =  coc-coL+  cog  m  -cbp<  0  such  that 
the  dissipative  terms  [^-contribution]  vanish,  Tm^n  (t)  contains  only 
dispersive  contributions, 

rm,„  (0  =  i  V  ^7, [S17] 

Xp  4A 

and  7m?n  is  defined  as  follows: 

Jm,n  =  Y  7 - — -  [S18] 

keBZ  (^L  -CDk-CDg^m+CDp) 


where  the  spin-exchange  coefficient  Jm/l  can  be  fully  controlled 
by  adjusting  Qa  and  cba.  We  have  the  following: 


~ 


Jm/l  - 


a,p 


-  0)g,n  +  d)a  —  Q)fj 


[S24] 


where  5(0)  =  1  and  8(x^  0)  =  0.  Note  that,  to  fully  control  Jm/l  at 
each  distance  rm/l,  we  need  to  introduce  enough  sidebands  to 
cover  all  of  the  energy  differences  cog^m  -  cog^n. 

If  the  characteristic  energy  scale  of  the  spin  Hamiltonian 
Hx y5o  is  much  smaller  than  the  minimum  energy  detuning 
8m  =  min{  | cog^n  -cog^n  \ }  between  different  sites,  that  is,  8w  ^ 
^2a^/(4A2),  the  time-dependent  processes  will  be  highly 


Hung  et  al.  www.pnas.org/cgi/content^hort/1 603^74^1  BUT|0N  A:  Distribution  approved  for  public  release. 


2  of  7 


off-resonant,  yielding  the  ideal  Hamiltonian  Hx y(t)  kH^ 0.  In 
practical  situations,  8m  will  be  a  limited  resource  because  of 
the  requirement  for  large  field  gradient  over  the  distance  of  a 
PCW  unit  cell.  Thus,  in  SI  Appendix  F:  Error  Reduction  and 
Analysis ,  we  discuss  errors  created  by  the  time-dependent 
processes,  and  strategies  to  minimize  the  errors. 

As  a  last  remark,  if  we  explicitly  write  down  the  phase  de¬ 
pendence  of  the  Raman  pump  in  Eq.  S23,  it  follows  that 
Jm,n  =  \Jm,n |c*kL'rw/z  (and  now  Jm,n  =Jn,m  for  kL-rm?n=0).  If  the 
illumination  is  not  perfectly  transverse,  that  is,  •  rm/l  ±  0, 
then  //xy,o  acquires  spatial-dependent,  complex  coupling 
coefficients: 


Ry(nj 2),  such  that  Ry(jr/2)(jxI^(jr/2)  =  oz,  to  transform  arbi¬ 
trary  XX  terms  into  desired  ZZ  interactions.  Spin  rotation  can 
be  realized,  for  example,  with  a  collective  microwave  driving 
ifmw  =  mw/2)oJL  +  h.c. ),  in  which  a  ^/2-microwave  pulse 

rotates  the  basis  { |g)„ ,  |s)„}  ->•  {(|g)„  +  \s)n)/A,  (-|g)„  +  \s)n)/\/2}. 

Thus,  an  Hxxz  Hamiltonian  can  be  simulated  using  the  following 
stroboscopic  evolution:  {Hxy , Hzz , Hxy , Hzz ,  •  •  •}  in  Nt  step  as 
schematically  depicted  in  Fig.  4.  The  unitary  evolution  in  each 
step  8t  =  t/Nt  is  given  by  the  following: 

e-i(HxY+Hzz )st  a e-iHxYdte-iHzzSt  _ l[HxY,Hzz}St2\  ,  jgjg] 


Hv, o=  t  +  [S25] 


where  we  see  that  the  leading  error  is  0(8t2).  When  repeating  this 
step  Nt  times,  the  leading  error  in  the  evolution  can  be  bounded 
by  the  following: 


which  gives  us  the  possibility  to  engineer  geometrical  phases  and 
nontrivial  topological  spin  models.  We  also  note  that,  in  the  above 
simple  expression,  we  are  assuming  ka  =  kL  being  a  constant  for 
all  different  sidebands  £la. 

Independent  Control  of  XX  and  YY  Interactions.  So  far,  we  have  been 
able  to  engineer  full  control  of  spin-exchange  or  XY  Hamilto¬ 
nians.  In  this  section,  we  will  show  how  by  slight  modification  of 
the  atomic-level  structure,  we  can  engineer  the  XX  and  YY  terms 
independently.  In  particular,  we  use  a  butterfly-like  structure  as 
depicted  in  Fig.  3 B,  where  there  are  two  transitions  coupled  to 
the  GMs,  that  is,  | g)  ++  \e)  and  | s)  \e),  and  two  different  mul¬ 
tifrequency  Raman  fields,  £lg(t)  and  Qs(t).  Assuming  that  we 
have  copropagating  beams  or  perfectly  transverse  illumination, 
that  is,  kL  •  rm,n  =  0,  we  can  adiabatically  eliminate  the  excited 
states  \e)  and  \e)  following  a  similar  procedure  as  in  the  pre¬ 
vious  section  (Eq.  Sll)  to  obtain  an  effective  light-matter 
Hamiltonian: 

//effjmM = -  ^  + h.c., 

k,n 

[526] 

where  we  assumed  that  £25(7) /As=  £lg(t)el(f,& /Ag  =  Q(t) / A,  Asg  = 
coe  -  a>L,sg,  and  (pgs  is  the  relative  phase  between  the  two  multi¬ 
frequency  Raman  fields  that  can  be  adjusted  at  will. 

Adiabatically  eliminating  the  photonic  modes  under  all  of  the 
approximations  that  we  used  in  the  previous  section,  we  arrive  at 
an  effective  Hamiltonian: 

Hxx,yy,0  =  \^m’n  (fgs  +  )  {°>lg  +  e  +  h.C.  j  , 

m,n>m  7 

[527] 

which,  depending  on  the  phase  (pgs,  can  drive  either  lor  7 
component,  that  is,  (o£s  ±  o™),  or  more  exotic  combinations  for 
general  (pgs.  Moreover,  if  the  two  pump  beams  are  not  copropa¬ 
gating,  they  create  spatial-dependent  phases  (pgs .  This  can  create 
site-dependent  XX,  YY,  or  XY  terms. 

Stroboscopic  Engineering  of  ZZ  Interactions.  One  way  to  engi¬ 
neer  Hzz  from  pure  XY  terms  is  to  notice  that  Rx(±ji/2) 
(o"<y?+o?o?)RU±n/2)=<%o?+onzo?  and  Ry(±n/2)(o"o™  + 
a^a^)R^,(±7i/2)  =a^a^  ±a^a^,  where  we  have  used  the  fol- 
lowing  notation  to  characterize  rotations  along  the  n  axis, 
that  is,  Rn(O)=ela'rv0/2.  This  is  particularly  useful  when  both 
XY  and  ZZ  interactions  have  the  same  coupling  strengths 
(e.g.,  the  Haldane-Shastry  model).  For  engineering  ZZ  terms 
in  a  generic  spin  Hamiltonian,  one  can  apply  spin  rotations 


<  \\[HxY,Hzz]\\t2 


[S29] 


for  Nt  1.  It  can  be  shown  (104)  that  higher-order  error  terms 
give  smaller  error  bounds.  Because  of  long-range  interactions, 
the  commutator  [Hxy,Hzz\  contains  up  to  N(N  -1)(N  -2) 
terms  that  are  different  from  0.1  Thus,  the  scaling  of  the  error 
in  the  limit  of  N  1  is  approximately  given  by  the  following: 


N{RJtf 
2-  Nt 


[S30] 


where  J  =  max [Jm,n\  is  the  largest  energy  scale  of  the  Hamiltonian 
we  want  to  simulate,  and  R  is  the  approximate  number  of  atoms 
coupled  through  the  interaction.  For  example,  if  Jm,n  is  a 
nearest  neighbor  interaction,  R  =  l.  If  Jm^  cx  l/\m  -n^,  then 
1/1^  T,  which  typically  grows  much  slower  than  N.  Be¬ 
cause  E2  cx  1  /Nt,  the  Trotter  error  in  Nt  steps  can  in  principle  be 
decreased  to  a  given  accuracy  s  by  using  enough  steps,  that  is, 
Nt  >N(RJt)2 /s. 

More  complicated  stroboscopic  evolutions  may  lead  to  a  more 
favorable  error  scaling  (64-66),  although  in  real  experiments 
there  will  be  a  trade-off  between  minimizing  the  Trotter  error 
and  the  fidelity  of  the  individual  operations  to  achieve  Hxy  and 
Hzz-  As  this  will  depend  on  the  particular  experimental  setup,  we 
will  leave  such  analysis  out  of  current  discussions.  For  illustra¬ 
tion,  we  will  only  consider  the  simplest  kind  of  stroboscopic 
evolution  that  we  depicted  in  Fig.  4. 

SI  Appendix  B:  Proper  Choice  of  Ground-State  Energy  Shifts 
in  2D  Models 

For  generic  2D  lattice  models,  we  need  to  introduce  ground-state 
energy  shifts  between  sites  (m,n)  separated  by  rm/l  =  rm  -rn  to 
engineer  the  interaction  between  them.  The  energy  shifts  need  to 
be  unique  for  specific  site  separation  vector  rm  n,  but  should  be 
independent  of  rm  to  preserve  translational  invariance,  which  is 
generally  required  for  spin  lattice  models.  To  do  this,  we  can 
introduce  a  linear  magnetic  field  gradient  VB  •  a/  =  Bi,  where  a / 
(/  =  1,2)  are  the  Bravais  vectors  of  a  unit  cell.  We  require  that  the 
ratio  q=B\/B2  of  B-field  gradients  along  two  Bravais  vectors 
be  an  irrational  number  such  that,  for  any  rn  =  (nx,ny)  with 
nx,nyG Z,  cog^=fiB{nxxq  +  ny)B2  is  a  unique  number.  As  a  re¬ 
sult,  each  sideband  can  only  induce  a  resonant  interaction  at  a 
specific  separation  rm  n.  Moreover,  we  also  need  to  ensure  that 
there  exists  no  rm  such  that  \cog/n  -cog,n \  <  Jm,n  that  can  lead  to 


+More  precisely,  [HXY,Hzz}  =  -2YZ^Zn^.n^gfnsg-  where  &m,n  =  Eq*m.n(Jq,m-Jq.n)°l 
Assuming  a  highly  frustrated  state,  we  simply  estimate  that  |(A^n)j  <  7/?  and 
\\[Hxy,Hzz]\\<2{JR)2. 
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significant  time-dependent  terms  in  Eq.  6  of  the  main  text.  In 
general,  for  a  finite-size  system,  this  situation  can  be  avoided. 

SI  Appendix  C:  Pump  Field  Configurations  for  Engineering  a 
Chiral-Flux  Square  Lattice  Model 

In  this  section,  we  explicitly  write  down  the  time-dependent  electric 
field  E(r n,t)  that  generates  the  desired  Raman  field  Q(t)  = 
(s|d  •  E*\e)  for  the  chiral-flux  Hamiltonian,  Eq.  15,  discussed  in  the 
main  text.  We  have  E(r n,t)  =z[sx{t)eM'rn  4-  sy{t)elky'rn]e~l(OLt,  where 
(s|d|e)  is  the  transition  dipole  moment.  For  the  field  propagating 
along  y,  the  amplitude  reads  as  follows: 

By 0)  =  eo  +  ei  (, e~isj  -  e ~is>‘)  +  e3  (e~^‘  +  e~w>‘) .  [S31] 

For  the  field  propagating  along  x,  we  similarly  require  the 
following: 

ex{t)  =  iCs  i  (e~is*‘  +  e~is>‘ )  +  e2  (e"'v  -  e~is^‘) .  [S32] 

Each  term  in  Eqs.  S31  and  S32  contributes  to  specific  sideband  in 
the  Raman  field  and  |£2a|  =  (s\d  -z\e)sa.  Pairing  individual  side¬ 
bands  to  the  main  Raman  field  introduced  by  the  leading  term 
eo  >  ^1,2,3  leads  to  the  desired  spin-exchange  interactions  as 
discussed  in  the  main  text. 

We  note  that  there  are  more  ways  other  than  Eqs.  S31  and  S32 
to  engineer  the  spin  Hamiltonian.  It  is  also  possible  to  introduce 
both  blue-detuned  (<5a>0)  and  red-detuned  ( 8-a  =  -8a )  side¬ 
bands  in  the  Raman  field  to  control  the  same  spin-exchange 
term.  That  is, 

Jm,n  —  J  j^o(rn)-^cr(rm)  T Xq  (rm)^L_a(rn)J  ,  [S33] 

which  has  contributions  from  Xa  and  X_a  of  blue  and  red  side¬ 
bands,  respectively.  Arranging  both  sidebands  with  equal  ampli¬ 
tudes  can  lead  to  equal  contributions  in  the  engineered  coupling 
coefficient.  This  corresponds  to  replacing  frequency  shifts  e~l5at 
in  Eqs.  S31  and  S32  with  amplitude  modulations  cos  5 at.  We  may 
replace  the  fields  by  the  following: 


(\Q.\/4A)[e~myJl  +  f £-*(%+!>],  formed  by  two  fields  propagat¬ 
ing  along  y  and  x,  respectively.  If  both  fields  have  the  same 
amplitude  (f  =  1),  they  either  add  up  or  cancel  completely, 
depending  on  whether  nx  +  ny  is  odd  or  even.  The  resulting 
coupling  rate  is  real  with  amplitude  as  follows: 

Jmfl  =JX ox;  =  |  [1  -  (-1)"^]  ,  [S37] 

and  vanishes  exactly  in  a  checkerboard  pattern.  If  one  applies  the 
same  trick  toward  NN  coupling  along  y,  but  with  f  ^  1,  the  cou¬ 
pling  amplitude  also  modulates  in  a  checkerboard  pattern.  Essen¬ 
tially,  all  three  NN  terms  around  a  brick  wall  vertex  can  be 
independently  controlled,  opening  up  further  possibility  to  engi¬ 
neer,  for  example,  Kitaev’s  honeycomb  lattice  model  (87,  88). 

•  Complex  NNN  and  NNNN  couplings:  The  NNNN  terms  can 
be  similarly  generated  by  using  a  sideband  of  detuning  2dy , 
andX2_y(r„  +  Ar2};)  =  ^[e~l(ny+T}n  +  i£e~mx!r],  formed  by  two  ini¬ 
tially  n/2  out-of-phase  fields  propagating,  respectively,  along  y 
and  x.  The  same  trick  can  be  used  for  NNN  couplings  using 
sidebands  of  detunings  8^,8^*  respectively.  The  coupling 
phase  (p  can  be  arbitrarily  controlled  by  the  amplitude  ratio 
f ,  as  we  showed  in  Eq.  16  in  the  main  text. 

To  engineer  this  model,  again  only  two  pump  beams  can  in¬ 
troduce  all  components  required  in  the  Raman  field.  We  can 
write  down  the  x-,  y-propagating  fields  with  amplitude  modula¬ 
tions  (that  is,  with  equal  blue-  and  red-sideband  contributions); 
for  the  field  propagating  along  y,  we  have  the  following: 


Sy(t)  =  so  +  s\  ^  cos  8xt  -  cos 

+  82  (cos  28yt  -  COS  8xyt  -  COS  8xy*  t) , 


[S38] 


and,  for  the  field  propagating  along  x,  we  require  the  following: 
ex(t)  =  y cos  8xt  -  i£s2  (cos  28yt  +  cos  8 -  cos  8^4) .  [S39] 


sy  (t)  =  s0  +  y  (cos  8xt  -  cos  8yt)  +  y  (cos  2 8xt  +  cos  2 8yt) ,  [S34] 

ex(t)  =  iC  (cos  &  +  cos  &yt)  +  y  ^cos  8xyt  -  cos  SxyJ^j .  [S35] 

SI  Appendix  D:  Pump  Field  Configurations  for  Engineering  a 
Topological  Spin  Model  in  a  Brick  Wall  Lattice 

In  this  section,  we  describe  in  detail  how  to  engineer  Haldane’s 
topological  spin  model,  Eq.  17  of  the  main  text, 

H=hJ2  (*»*»  +  h-c)  +  ^2  E  {eiKnol°n  +  h.c.) ,  [S36] 

(, m,n )  {m,n} 

in  a  brick  wall  configuration.  Introducing  a  strong  pump  field  of 
amplitude  £2o>  propagating  along  y,  we  can  engineer  the  interac¬ 
tion  terms  one-by-one  as  the  following: 

•  Uniform  NN  coupling  along  y:  This  term  can  be  realized 
with  a  sideband  of  detuning  8y  and  Xy(rn  +  Ary)  =  -(|£2i|/ 
2A )^-?(%+1)^  pairing  with  the  strong  pump  field  Xo(r„)  = 
(|£2o|/2A )e~lnyn,  which  propagates  along  y.  The  resulting  cou¬ 
pling  coefficient  is  ti=J\Xo\\Xy\.  This  term  can  also  be  ex¬ 
tended  to  engineer  nonuniform  coupling  coefficients;  see 
the  following  discussion. 

•  Checkerboard-like  NN  coupling  along  x:  We  introduce 
a  sideband  of  detuning  8X  and  amplitude  Xx(rn  +  Arx)  = 


One  may  similarly  replace  amplitude  modulations  cos  8at  by  fre¬ 
quency  modulation  e~l5at  to  engineer  the  spin  model,  as  in  Eqs. 
S31  and  S32. 

SI  Appendix  E:  PCW  and  Pump  Field  Configurations  for 
Engineering  an  XXZ  Spin  Hamiltonian  with  Mr"1  Interaction 

We  describe  how  to  engineer  an  XXZ  Hamiltonian,  which  is 
typically  written  as  follows: 

Hxxz  =  -B  E  <+E  |r  ZZ  1 1  [COS(6,)'^C4 

n  n<m\rn  Im  I 

+  sin  (0)  (<<4 +*£<)], 

[S40] 

where  an  effective  magnetic  field  B  that  controls  the  number  of 
excitations  can  be  introduced  by  the  detuning  of  addressing 
beams;  the  parameter  6  determines  the  relative  strength  between 
the  ZZ  and  XY  interactions. 

We  can  simulate  any  rj  and  6  in  a  2D  PCW  as  follows: 

•  First,  we  require  the  decay  length  scale  £  to  be  sufficiently 
large  such  that  GM  photon  coupling  strength  7(rm?n)  is  only 
limited  by  the  energy  spread  in  the  low-dimensional  reservoir. 
For  example,  in  2D  with  quadratic  dispersion  as  depicted  in 
Fig.  ID  of  the  main  text,  one  can  simulate  any  interaction  that 
decays  faster  thanXo(|rm,„|/£)  «log(f/|rmjn|),  where  Kq(x)  is  the 


Hung  et  al.  www.pnas.org/cgi/content^hort/1 603^74^1  BUT|0N  A:  Distribution  approved  for  public  release. 


4  of  7 


modified  Bessel  function  of  the  second  kind  (36,  37)  (SI  Appendix 
A:  Complete  Derivation  of  Final  Time-Dependent  Hamiltonian). 

•  Then,  to  simulate  the  |r„  dependence,  we  introduce  a 

linear  magnetic  field  gradient  WB  or  linear  ground-state  en¬ 
ergy  shifts  as  described  in  SI  Appendix  B:  Proper  Choice  of 
Ground- State  Energy  Shifts  in  2D  Models. 

•  Thus,  as  we  did  in  the  previous  discussions,  we  introduce  a 
strong  pump  field  of  amplitude  £2o  together  with  Nd  auxiliary 
fields  Qa  of  detuning  8a  to  cover  all  of  the  different  separa¬ 
tions  rmp.  For  example,  to  simulate  a  square  lattice  of  nsxns 
(=N)  atomic  spins,  we  can  see  that  the  number  of  different 
distances  grows  as  Nd  =  (ns(ns  +  1)  -2)/2,  linearly  propor¬ 
tional  to  the  number  of  atoms  Nd  oc  A. 

•  Finally,  the  parameter  6  that  determines  the  ratio  between  ZZ 
and  XY  interaction  can  be  controlled  by  using  different  pump 
intensities  when  doing  the  stroboscopic  evolution. 

SI  Appendix  F:  Error  Reduction  and  Analysis 
Engineering  a  Two-Step  Hamiltonian.  As  derived  in  SI  Appendix  A: 
Complete  Derivation  of  Final  Time-Dependent  Hamiltonian , 
the  resulting  time-dependent  couplings  introduced  by  a  mul¬ 
tifrequency  pump  is  H(t)  =  ^2pHpeip5t,  where  Hp  represents  the 
part  in  the  Hamiltonian  that  oscillates  with  frequency  p8.  The 
time-dependent  Hamiltonian  has  a  period  T  —  2n/8,  and  it  can  be 
shown  that  the  effective  Hamiltonian  that  repeats  every  time  period 
T\  that  is,  e~iHcffT,  where  the  measurements  would  take  place  in  an 
experiment,  is  given  by  the  following  (99): 


#eff,l 


[HP,H-P\ 

P 


1  ^\[Hp,H0],H_p\  +  [[H_p,H0],Hp] 
282  ^  p2 


[S41] 


This  means  that,  if  we  apply  the  multifrequency  pumps  just  as 
we  explained  in  the  previous  sections,  the  leading  error  would 
be  on  the  order  of  J2  /  8,  where  J  is  the  interaction  strength  that 
we  want  to  simulate.  However,  we  also  observe  that  the  leading 
error  term  '  [Hp,H-p]/(p8)  vanishes  if  Hp  =  ±H_p.  In  other 
words,  first-order  error  vanishes  if  Hp  is  either  symmetric  or  anti¬ 
symmetric  under  time  reversal  operation  2Y.  Although  the  original 
Hamiltonian  H(t)  does  not  necessarily  possess  such  symmetry,  it 
is  possible  to  introduce  a  two-step  periodic  operation  F^step  = 
{H,  2TH,H,2TH,  . . .  }  to  cancel  the  first-order  error  while  keep¬ 
ing  the  time-independent  part  F/2step,o  =Ho  identical.  This  reduces 
the  leading  error  to  the  order  of  J3  /82. 

To  achieve  the  time  reversal  operation,  we  must  reverse  the 
phase  of  the  driving  lasers,  as  well  as  the  sign  of  the  energy  offsets 
between  the  atoms.  Specifically,  we  can  engineer  the  two-step 
Hamiltonian  as  follows: 

•  After  applying  a  proper  magnetic  (or  Stark  shift)  gradient  to 
ensure  a  position-dependent  energy  shift,  we  apply  sidebands 
using  either  frequency  or  amplitude  modulations  to  engineer 
the  Hamiltonian  Hq  following  Eq.  S23. 

•  After  a  time  T  =  2n/ 8,  we  flip  the  sign  of  the  gradient.  In  ID,  if 
cogp  -  cog^m  =  \n-m\8,  then  we  switch  to  o)gp  -  ojgm  =  —\n—m\8. 
Meanwhile,  we  also  reverse  the  propagation  direction  of 
the  Raman  fields  such  that  Xa  -+X*.  As  a  result,  all  of  the 
time-dependent  Hamiltonians  Hp,  \/p  ±  0,  become  H_p  in  the 
second  step;  whereas  the  time-independent  Hamiltonian 
Ho  remains  identical.  After  holding  for  an  evolution  time 
T  =  2n/8,  we  repeat  the  first  step. 

To  formally  prove  that  the  above  idea  leads  to  smaller  error,  we 
can  write  our  two-step  periodic  Hamiltonian  as  follows: 


FF2step (0  =H(t)T(t)  +H(—t)T(t  -  T),  [S42] 


where  T(t)  is  a  periodic  square-wave  envelope,  controlling  the  on 
and  off  of  H(t)  at  time  interval  [0,  T\  within  a  period  2 T.  T(t )  can 
be  expanded  as  follows: 


m= 


1  i 

~  +  — 

2  m 


E 

m  odd 


(e" 


;mdt/ 2  _ 


0-imdt, 


72) 


m 


[S43] 


where  m  =  1,3,  . . ..  Plugging  the  expansion  H(t)  =  "f2pHpeipdt  into 
Eq.  S42,  we  now  have  the  following: 


FF2step  (0  —  rs 


V*  1+w  E 


^imdt/2  _  g-im8t/2^j 


E 


m  odd 

gimdt/ 2  _  g-imSt/ 2) 


[S44] 


Writing //2step(0  =  EA  eip(5/2)t^  tpe  pourier  components  Hp  of 
the  two-step  periodic  Hamiltonian  are  as  follows: 


Hp(  even)  —  ^  (Hp/2  +h_p/2), 


[S45] 


Hp(odil)  -  ^  E  m  (H(p-m)/ 2  H-(p-m)/2  +H(p+m)/2  H-{p+m)/ 2)- 

m  odd 

[S46] 


So  we  have  Hp  =  (-VfH_p  and  the  leading  order  error  in  Eq.  S41 
vanishes.  According  to  the  Floquet  theory,  we  arrive  at  an  effec¬ 
tive  time-independent  Hamiltonian  FFeff,2  at  every  time  interval 
T2  =  An/ 8, 


//eft, 2  =Ho  +He rr,2  »//0  +4  E  (-1  fl^EllFAA.  [S47] 

8  p  P 

Time-Dependent  Stark  Shifts  in  the  Error  Analysis.  In  the  previous 
discussions,  we  have  dropped  the  contribution  of  the  time- 
dependent  Stark  shifts: 


mp—l 

//acW=-EESR 

n  a>p 


QMl 


2A 


<£> 


[S48] 


where  ooa^  =  oba~  a>/3-  In  the  following,  we  discuss  its  role  in  the 
effective  Hamiltonian,  using  the  Floquet  error  analysis. *  The 
Fourier  coefficients  of  the  Stark  shifts  can  be  written  as  follows: 

H*c,P=J2AP^’  [S49] 

n 


where  the  on-site  amplitude 

a;  =  -A  Y^x<x;s{&„,p  -Ps) ;  [S50] 

Ap  may  be  site  dependent  if  the  phase  differences  between  the 
Raman  fields  Qa  vary  across  sites.  Comparing  Eq.  S50  with  Eq.  8 
of  the  main  text,  we  see  that  A"  ~0(JA/J)  may  be  even  larger 


*Here,  we  discuss  a  special  case  where  only  |s)  state  is  shifted.  For  a  general  case  when 
both  |g)  and  |s)  states  are  shifted,  see  later  discussion. 
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than  the  engineered  interaction J  =  ma x\Jm/l]  if  A  >  7.  In  the  two- 
step  driving  scheme,  we  replace  the  amplitude  A ”  by  A  according 
to  Eqs.  S45  and  S46. 

As  discussed  in  the  previous  section,  with  two-step  driving, 
the  leading  Stark  shift  error  contribution  only  appears  in  the 
second  order: 


Time-Dependent  Stark  Shifts  in  the  Butterfly  Scheme.  We  discuss  the 
error  contribution  due  to  time-dependent  Stark  shifts  in  the 
butterfly  scheme  for  independent  control  of  XX  or  YY  interac¬ 
tions  (as  well  as  ZZ  interaction  in  stroboscopic  evolution).  The 
time-dependent  Stark  shift  is  as  follows: 


Ht 


err, 2 


4  ,  y,  [  \Hp  +  Hacp,  Hq\  ,  Hp  +  Hac,p\ 

<52V  T  p 2 


[S51] 


where  Hacp  are  the  Fourier  coefficients  of  the  two-step  Stark- 
shift  Hamiltonian.  Only  the  following  nested  commutators, 
[[Ha.c,p>H()\>Hac,p\)  [[Hac,p?Ho\?-Hp\>  and  [[Hp,Ho\,Hacp\  are  re¬ 
lated  to  the  time-dependent  Stark  shifts,  and  should  be  evaluated 
in  various  configurations  as  follows. 

Generic  Hamiltonians  with  translational  invariance.  By  translational 
invariance,  we  mean  that  there  are  no  site-dependent  spin  inter¬ 
actions,  and  the  spin-exchange  coefficients  remain  identical  as  we 
offset  the  spin  index  by  one  or  more.  This  means  that  all  components 
in  the  pump  field  should  drive  the  system  with  uniform  optical  phases 
as  in  the  Haldane-Shastry  model  discussed  above.  The  resulting 
Fourier  coefficients  Hdcp  =Ap^jnonss  would  have  identical  effect  on 
all  spins  (Ap  being  a  constant  amplitude).  As  a  result,  the  above- 
mentioned  commutators  vanish,  suggesting  that  the  error  by  Hac  (t) 
averages  out  to  zero  in  the  Floquet  picture.  In  the  butterfly  scheme, 
however,  both  |g)  and  | s)  states  are  pumped  and  they  may  be  shifted 
differently.  This  leads  to  slight  modifications  in  the  engineered  XX 
and  YY  terms  (see  later  discussion). 

Models  containing  sublattices.  For  topological  models  that  contain 
sublattices,  as  in  our  examples,  the  pump  fields  are  not  perfectly 
transverse  and  Stark  shifts  are  site  dependent,  resulting  in 
nonvanishing  error.  However,  we  also  note  that  the  coupling  co¬ 
efficients  appearing  in  the  commutators  [[Hacp,Ho\,Hacp\  and 
[{Hp,H0],HK,p\({{HaCj:,Ho\,Hp])  are  on  the  order  of/|£2|4/A2  and 
/2|£2|2/A,  respectively.  For  realistic  PCW  realizations,  we  may  set 
3/ A  >  0(1).  Because  Ap/Jm/i  ~0(A/J),  the  energy  scales  of  the 
commutators  are  all  on  the  order  </3,  and  the  energy  scale  in 
He n-,2  will  be  <J3/S2.  Therefore,  for  <5  »/,  the  Stark  shift  terms 
may  be  ignored. 

Stark  shift-dominated  regime.  It  may  be  possible  that  our  sublattice 
models  be  purposely  driven  with  large-amplitude  pumps  such  that 
|£2|2/A  >  5.  Stark  shift  contributions  would  become  important  in 
the  resulting  spin  dynamics.  However,  if  we  choose  a  large  pump 
detuning  A  >/,  the  dominant  error  contribution  [recall  that 
/~(//A)(|£2|2/A)]  comes  from  the  [[Hacp,H o\,Hacp\  term,  which 
can  be  written  in  the  following  simple  form: 


Hen, 2  A  £^  [[H^HolH^j 
0  p  y 

«  £^m,n  (jmsiQ^otg  +  ft-c)  ■ 

m,n 


[S52] 


Flere,  the  site-dependent  Fourier  coefficient  is  defined  as  HXJ)  = 
Y.M  and^m,n  =  Ep((4(-ir1)/52F2)M;-^)V0in  gen- 
eral  for  interactions  Jm/l  that  are  not  translationally  invariant.  In 
a  special  case  that  only  two  sublattices  are  present,  as  in  our 
examples,  we  note  that  Am^n  may  only  depend  on  the  distance 
rm/l  and  is  site  independent.  This  “error”  term  would  then  uni¬ 
formly  modify  the  XY  coupling  strengths  to  a  modified  value: 


Jm,n  =  (1  +A  m,n)Jm,n • 


[S53] 


The  next  leading  order  errors  are  due  to  [[Hp^H^.H^^/d2  and 
[[Hacp,Ho\,Hp\/82  terms,  which  are  on  the  order  of  J2\tl\2  /  (A52) 
and  are  a  factor  of  ~J/ A  smaller  than  the  leading  Stark  shift 
contribution.  This  suggests  we  can  always  increase  the  detuning 
A,  while  keeping  |Q|/A  constant,  to  reduce  the  error  contribution. 
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where  cbaj  =  doa-  (bp.  The  Fourier  coefficients  of  the  Stark  shifts 
can  be  written  as  follows: 


^  A  V 


[S55] 


where  A"  =  -  ^(KQ.aQ*JAA)8(d)a£-pd),  and  we  have  used  that 

a+B 

fact  that  |£\;a|/Ag  =  |^a|/A5  =  |£^a|/A.  In  two-step  driving,  we 

replaced"  with  Ap  according  to  Eqs.  S45  and  S46.  We  note  that 
states  |g)  and  | s)  can  be  driven  differently  when  Ag  ±  As.  This 
leads  to  different  error  comparing  to  simple  Raman  driving  only 
on  one  of  the  state. 

In  two-step  driving,  the  leading  Stark  shift  error  contribution 
appears  in  the  second  order: 


H&xr,2 


4  ^  y?  \  \Hp  o]  ,Hp  -\-Hac,p\ 


[S56] 


where  Hacp  are  the  Fourier  coefficients  of  the  two-step  Stark 
shift  Hamiltonian.  To  simplify  the  discussion,  we  discuss  XX  in¬ 
teraction  with  <pgs  =  0  and  the  following: 


H<>  =Hxx.<)  =  £  /«*(<£  +<g)(<g  +  oj). 


[S57] 


For  illustration,  we  evaluate  the  following  nested  commutators, 
[ [33 dCp ,  Hq\ ,f3dcp\,  to  access  the  error  contribution  due  to  time- 
dependent  Stark  shifts.  We  find  the  following: 


P  ac  ,p?  H0],H 

ac  ,p 


N  I 

EJm,n  /  .  . 

Tr(A*-A* 

m,n>m  ^ 

+°Kr)  +  (pP  -ap 


)2[(4m+^”)2(« 

)'(«+«) 


[S58] 
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We  therefore  see  that,  as  long  as  Ag  =  As,  there  is  no  error  due 
to  time-dependent  Stark  shifts  as  a  result  that  both  | g)  and  | s) 
states  are  shifted  exactly  the  same  way.  However,  one  may  find 
that  the  atomic-level  structure  for  a  butterffy  scheme  dictates 
that  Ag^As.  As  a  result,  we  find  [[H^,Ho],Hacj)]  ^0  even  for 
translational-invariant  models,  where  A  =Ap  are  independent  of 
sites.  This  is  in  contrast  to  the  case  of  XY  Hamiltonians  with  the 
Raman  field  driving  only  one  ground  state.  It  is,  however,  still 
possible  to  minimize  the  error  contribution  by  driving  with  J>  A 
and  5  >  /,  as  the  criterion  needed  for  sublattice  models  (see 
main  text). 
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In  the  case  of  strong  driving  with  Ag^As  and  A  >/,  the  term  is  in  fact  also  driving  XX  and  YY  interactions.  It  is  therefore 

commutator  [[H^p,Hq\,H^p]  can  be  made  as  the  dominant  desirable  to  take  the  Stark  shift  contribution  into  account  while 

error  contribution.  Interestingly,  as  seen  from  Eq.  S59,  this  error  finding  the  proper  pump  fields  to  create  the  target  Hamiltonians. 
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Tailoring  the  interactions  between  quantum  emitters  and  single 
photons  constitutes  one  of  the  cornerstones  of  quantum  optics. 
Coupling  a  quantum  emitter  to  the  band  edge  of  a  photonic  crystal 
waveguide  (PCW)  provides  a  unique  platform  for  tuning  these 
interactions.  In  particular,  the  cross-over  from  propagating  fields 
£(x)  oc  e±lkxX  outside  the  bandgap  to  localized  fields  E{x)oce~Kx^ 
within  the  bandgap  should  be  accompanied  by  a  transition  from 
largely  dissipative  atom-atom  interactions  to  a  regime  where  dis¬ 
persive  atom-atom  interactions  are  dominant.  Here,  we  experi¬ 
mentally  observe  this  transition  by  shifting  the  band  edge 
frequency  of  the  PCW  relative  to  the  D<|  line  of  atomic  cesium 
for  N  =  3.0  ±0.5  atoms  trapped  along  the  PCW.  Our  results  are 
the  initial  demonstration  of  this  paradigm  for  coherent  atom- 
atom  interactions  with  low  dissipation  into  the  guided  mode. 

quantum  optics  |  nanophotonics  |  atomic  physics 

Recent  years  have  witnessed  a  spark  of  interest  in  combining 
atoms  and  other  quantum  emitters  with  photonic  nano¬ 
structures  (1).  Many  efforts  have  focused  on  enhancing  emission 
into  preferred  electromagnetic  modes  relative  to  vacuum  emis¬ 
sion,  thereby  establishing  efficient  quantum  matter-light  inter¬ 
faces  and  enabling  diverse  protocols  in  quantum  information 
processing  (2).  Photonic  structures  developed  for  this  purpose 
include  high-quality  cavities  (3-7),  dielectric  fibers  (8-13),  me¬ 
tallic  waveguides  (14-16),  and  superconducting  circuits  (17-19). 
Photonic  crystal  waveguides  (PCWs)  are  of  particular  interest, 
because  the  periodicity  of  the  dielectric  structure  drastically  modifies 
the  field  propagation,  yielding  a  set  of  Bloch  bands  for  the  guided 
modes  (GMs)  (20).  For  example,  recent  experiments  have  shown 
superradiant  atomic  emission  because  of  a  reduction  in  group  ve¬ 
locity  for  an  atomic  frequency  near  a  band  edge  of  a  PCW  (21). 

A  quite  different  paradigm  for  atom-light  interactions  in  pho¬ 
tonic  crystals  was  proposed  in  the  works  in  refs.  22-25  but  has  yet 
to  be  experimentally  explored.  In  particular,  when  an  atomic 
transition  frequency  is  situated  within  a  bandgap  of  a  PCW,  an 
atom  can  no  longer  emit  propagating  waves  into  GMs  of  the 
structure.  However,  an  evanescent  wave  surrounding  the  atoms 
can  still  form,  resulting  in  the  formation  of  atom-photon-bound 
states  (26,  27).  This  phenomenon  has  attracted  new  interest  re¬ 
cently  as  a  means  to  realize  dispersive  interactions  between  atoms 
without  dissipative  decay  into  GMs.  The  spatial  range  of  atom- 
atom  interactions  is  tunable  for  ID  and  2D  PCWs  and  set  by  the 
size  of  the  photonic  component  of  the  bound  state  (28,  29).  Many- 
body  physics  with  large  spin  exchange  energies  and  low  dissipation 
can  thereby  be  realized  in  a  generalization  of  cavity  quantum 
electrodynamics  (CQED)  arrays  (30,  31).  Fueled  by  such  perspec¬ 
tives,  there  have  been  recent  experimental  observations  with  atoms 
(21,  32,  33)  and  quantum  dots  (34,  35)  interacting  through  the  GMs 
of  PCWs,  albeit  in  frequency  regions  outside  the  bandgap,  where 
GMs  are  propagating  fields. 

In  this  manuscript,  we  report  the  observation  of  collective  dis¬ 
persive  shifts  of  the  atomic  resonance  around  the  band  edge  of  a 
photonic  crystal.  Thermal  tuning  allows  us  to  control  the  offset  of 


the  band  edge  frequency  (ebe)  of  the  PCW  relative  to  the  frequency 
vdi  of  the  Di  line  of  cesium  (Cs).  In  both  the  dispersive  do¬ 
main  [i.e.,  vdi  outside  the  bandgap  with  electric  field  E(x)  oce±lkxX] 
and  reactive  regime  [i.e.,  edi  inside  the  bandgap  with  E(x)  (xe~Kx^], 
we  record  transmission  spectra  for  atoms  trapped  along  the  PCW, 
as  illustrated  in  Fig.  L4. 

To  connect  the  features  of  the  measured  transmission  spectra 
to  underlying  atom-atom  radiative  interactions,  we  have  developed 
a  formalism  based  on  the  electromagnetic  Green’s  function.  The 
model  allows  us  to  infer  the  peak  single-atom  frequency  shift  of 
the  atomic  resonance  /id  (Are)  and  GM  decay  rate  riD(ABE)  as 
functions  of  detuning  ABE  =  edi  -  ebe  between  the  atomic  edi  and 
band  edge  ebe  frequencies.  From  the  observation  of  superradiant 
emission  outside  the  bandgap,  we  infer  the  average  number  of 
trapped  atoms  to  be  N  =  3.0  ±  0.5,  as  described  in  ref.  21  and  SI 
Text.  (SI  Text  has  thorough  descriptions  of  the  design  and  charac¬ 
terization  of  the  PCW,  how  to  obtain  the  attenuation  coefficient 
and  the  band  edge  position  of  the  PCW,  how  to  generate  the 
atomic  spectra  fits,  and  the  measurements  of  atomic  decay.)  For 
frequencies  inside  the  bandgap  (ABE  =  50  GHz),  the  ratio  of  dis¬ 
sipative  to  coherent  rates  is  H  =  Tyd/Jyd  =  0.05  ±0.17  because  of 
the  exponential  localization  of  the  atomic  radiation  in  the  bandgap. 
For  comparison,  the  prediction  for  our  system  from  CQED  models 
alone  is  T^cqed  =  0.30  ±  0.04.  Other  than  yielding  a  more  favorable 
ratio  between  coherent  and  dissipative  GM  rates,  PCWs  offer 
significant  advantages  compared  with  conventional  cavities  as 


Significance 

In  recent  years,  there  has  been  considerable  effort  to  bring 
ultracold  atoms  into  the  realm  of  nanophotonics.  Nanoscopic 
dielectric  devices  offer  unprecedented  opportunities  to  engi¬ 
neer  novel  capabilities  for  the  control  of  atom-photon  inter¬ 
actions.  In  particular,  photonic  crystals  are  periodic  dielectric 
structures  that  display  a  photonic  bandgap  where  light  cannot 
propagate  and  provide  a  new  setting  for  coherent  photon- 
mediated  interactions  between  atoms  with  tunable  range.  Here, 
we  report  the  initial  observation  of  cooperative  atom-atom  in¬ 
teractions  around  the  band  edge  of  a  photonic  crystal  wave¬ 
guide.  Our  experiment  opens  the  door  to  fascinating  scenarios, 
such  as  exploring  many-body  physics  with  large  spin  exchange 
energies  and  low  dissipation. 
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platforms  for  atom-light  interfaces.  First,  the  range  of  interaction 
in  a  PCW  is  tunable,  ranging  from  effectively  infinite  to  nearest 
neighbor  (28,  29,  36),  in  contrast  to  the  fixed  infinite  range  of  a 
cavity.  Second,  because  of  the  multimode  nature  of  PCWs,  one 
can  use  different  GMs  as  different  interaction  channels  to 
which  the  atoms  simultaneously  couple. 

Alligator  PCW 

Fig.  L4  provides  an  overview  of  our  experiment  with  atoms 
trapped  near  and  strongly  interacting  with  the  transverse- 
electric  (TE)  mode  of  an  alligator  PCW.  The  suspended  silicon 
nitride  (SiN)  structure  consists  of  iVceiis  =  150  nominally  identical 
unit  cells  of  lattice  constant  a  =  370  nm  and  is  terminated  by  30 
tapering  cells  on  each  side,  as  shown  in  the  SEM  images  in  Fig.  IB. 
The  tapers  mode-match  the  fields  of  the  PCW  to  the  fields  of 
uncorrugated  nanobeams  for  efficient  input  and  output  coupling. 
Design,  fabrication,  and  characterization  details  are  described  in 
refs.  21,  32,  and  33.  Fig.  1C  shows  the  nominal  cell  dispersion 
relations  for  the  TE  (polarized  mainly  along  y )  and  transverse- 
magnetic  (TM)  modes  (polarized  mainly  along  z).  After  release  of  the 
SiN  structure  from  the  silicon  (Si)  substrate,  a  low-power  CF4  etch  is 
used  to  align  the  lower/“dielectric”  TE  band  edge  (vbe)  to  the  Cs  Di 
transition  (vDi)-  The  TM  mode  has  band  edges  far  detuned  from  the 
both  the  Cs  Di  and  D2  lines.  In  our  experiment,  the  TE  mode  is  used 
to  probe  the  atoms,  whereas  the  TM  mode  with  approximately 
linear  dispersion  serves  to  calibrate  the  density  and  trap  properties. 

To  better  understand  atomic  interactions  with  the  PCW,  it  is 
helpful  to  visualize  the  spatial  profile  of  the  fields  generated  absent 
atoms,  when  light  is  input  from  one  end.  Fig.  24  shows  the  mea¬ 
sured  intensity  along  the  length  of  the  PCW  as  a  function  of  probe 
detuning  5be  =  vp-  vbe  around  the  band  edge,  where  vp  is  the 


probe  frequency.  The  intensity  was  measured  by  imagmg  weak 
scatterers  along  the  length  of  the  alligator  PCW  that,  after  cali¬ 
bration,  serve  as  local  probes  of  the  intensity  (SI  Text).  Fig.  2 B 
shows  the  corresponding  finite  difference  time  domain  (FDTD) 
simulated  intensity  (37).  In  both  images,  resonances  appear  at 
vp  =  vi?2,3  because  of  the  weak  cavity  formed  by  the  reflections  of 
the  tapers.  The  spatial  modulation  of  the  intensity  at  the  resonances 
caused  by  the  cavity  effect  is  approximated  by  \E(x)\2  « cos2 (5 kx  x), 
where  6kx  =  n/a  —  kx  is  the  effective  wavevector  near  the  band  edge. 
The  nth  resonance  at  frequency  vn  is  such  that  8kx  =  n%/L ,  where 
L  is  the  effective  length  of  the  PCW  (including  field  penetration 
into  the  tapers).  Fig.  2C  shows  a  plot  of  \E(x)\  for  a  probe  input 
at  frequency  vp  =  v\  at  the  first  resonance.  Inside  the  bandgap 
(Abe  >  0),  the  field  is  evanescent,  and  5 kx  —  Ikx.  Fig.  2D  plots  | E(x)  \ 
for  probe  frequency  vp  =  vbg  inside  the  bandgap  and  shows  the 
exponential  decay  of  the  intensity.  Using  a  model  for  the  field  in  a 
finite  photonic  crystal  (SI  Text),  we  fit  the  measured  intensity  for  each 
frequency  in  Fig.  2  A  and  B  and  extract  5 kx  and  kx,  thereby  obtaining 
the  dispersion  relations  shown  in  Fig.  2 E.  Importantly,  we  determine 
the  band  edge  frequency  for  the  actual  device  to  be  eBe  -v\  =  133  ±  9 
GHz  relative  to  the  readily  measured  first  resonance  at  v\,  which  is  in 
good  agreement  with  the  FDTD-simulated  result  of  135  GHz. 

Both  v\  and  vbg  are  relevant  to  our  measurements  of  trans¬ 
mission  spectra  with  trapped  atoms.  The  presence  of  a  “cavity” 
mode  at  v\  implies  that  the  emission  of  an  atom  with  transition 
frequency  vdi  =  v\  will  generate  a  field  inside  the  PCW  with  an 
analogous  spatial  profile  to  that  of  the  cavity  mode,  as  shown  in 
Fig.  2 C.  By  contrast,  atomic  emission  in  the  regime  with 
vdi  =  vBg  within  the  bandgap  will  excite  an  exponentially  local¬ 
ized  mode  centered  around  the  atomic  position  jca,  as  illustrated 
in  Fig.  IF. 


Fig.  1.  Description  of  the  alligator  PCW.  (>A)  Atoms  are  trapped  above  the  PCW  in  an  optical  dipole  trap  formed  by  the  reflection  of  a  near-normal 
incidence  external  beam  (21).  The  orange  cylinder  represents  the  confinement  of  the  atoms,  which  is  Axa  ~  ±6  ^im  along  the  axis  of  the  device  and 
A yA  ~  A zA  ~  ±30  nm  in  the  transverse  directions  (SI  Text).  The  three  green  spheres  represent  trapped  atoms  that  interact  radiatively  through  the  fun¬ 
damental  TE  GM,  polarized  mainly  along  y.  The  decay  rate  for  a  single  atom  into  the  PCW  is  r1D  (red  arrows),  and  the  decay  rate  into  all  other  modes  is  F 
(wavy  red  arrow).  ( B )  SEM  images  of  portions  of  the  tapering  and  PCW  sections.  The  suspended  SiN  device  (gray)  consists  of  1 50  cells  and  30  tapering  cells 
on  each  side.  The  lattice  constant  is  a  =  370  nm,  and  thickness  is  185  nm.  (C)  Calculated  band  structure  of  the  fundamental  TE  (solid)  and  TM  (translucent) 
modes  using  an  eigenmode  solver  (38)  and  the  measured  SEM  dimensions,  which  are  modified  within  their  uncertainty  to  match  the  measured  bands.  The 
black  curves  represent  the  Bloch  wavevector  kx  (lower  axis).  The  red  curves  show  the  attenuation  coefficient  kx  of  the  field  for  frequencies  in  the  bandgap 
(upper  axis)  and  are  calculated  by  means  of  an  analytical  model  (SI  Text).  The  dotted  lines  mark  the  frequencies  of  the  Cs  Di  (vDi  =335.1  THz)  and  D2 
(vD2  =  351.7  THz)  transitions.  The  dielectric  band  edge  is  indicated  as  vBe-  The  pink  shaded  area  represents  the  TE  bandgap.  The  gray  shaded  area  rep¬ 
resents  the  light  cone. 
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Fig.  2.  Characterization  of  the  alligator  PCW.  (/A)  Measured  and  ( B )  calculated  electric  field  magnitudes  along  the  PCW  as  functions  of  position  x  along 
the  PCW  and  probe  detuning  5BE=vp-vBE  relative  to  vBe  for  the  dielectric  band  edge.  (C  and  D)  GM  intensity  |£(x)|2  along  PCW  at  two  different  fre¬ 
quencies:  (C)  vi  for  the  first  cavity  resonance  showing  a  resonant  "supermode"  and  (D)  vBG  inside  the  bandgap  displaying  exponential  decay  (/VcensKxa  =  2.0 
at  vBG).  For  clarity,  the  number  of  cells  of  the  nominal  and  tapering  sections  is  decreased  by  a  factor  of  five,  and  the  Bloch  periodicity  (a  =  370  nm),  al¬ 
though  present,  is  not  shown  in  the  intensity.  The  orange  ovals  represent  the  confinement  of  the  atoms  in  the  optical  trap  above  the  PCW,  which  is 
Axa  ~  ±6  |am  along  the  x  axis  of  the  device  and  A yA  ~  ±30  nm,  with  a  PCW  gap  width  of  220  nm.  (£)  Dispersion  relation  for  the  projected  wavevector  kx 
and  attenuation  constant  kx  vs.  probe  detuning  5Be  deduced  for  the  PCW  obtained  by  fitting  the  data  in  A  to  a  model  of  the  device  (5/  Text).  The  shaded 
pink  area  represents  frequencies  inside  the  bandgap.  (E)  Plot  of  the  exponentially  localized  emission  e_2K*lx_x^  from  an  atom  (green  sphere)  at  position  xA 
with  transition  frequency  vDi  =  vBG  inside  the  bandgap. 


Experiment 

Cs  atoms  are  trapped  above  the  surface  of  the  alligator  PCW,  as 
shown  in  Fig.  L4,  using  a  similar  experimental  setup  to  that 
reported  in  ref.  21.  As  described  in  more  detail  in  ref.  21,  the  decay 
rate  T\D  into  the  GM  is  exponentially  sensitive  to  the  trap  position 
above  the  surface  of  the  alligator  PCW.  Our  calculations  and 
measurements  of  r m  agree  with  COMSOL  simulations  (38)  of  the 
trap  position,  and  thus,  we  are  able  to  determine  that  the  Cs  atoms 
are  trapped  145  ±  15  nm  above  the  surface  of  the  alligator  PCW. 
Atoms  are  cooled  and  trapped  in  a  magneto-optical  trap  (MOT) 
around  the  PCW  and  then  loaded  into  a  dipole  trap  formed  by  the 
reflection  from  the  device  of  a  frequency  red-detuned  side  illumi¬ 
nation  (SI)  beam.  The  SI  beam  has  a  waist  of  50  pm,  and  the  po¬ 
larization  is  aligned  along  the  x  axis  for  maximum  reflection  from 
the  PCW.  We  measure  a  \/e  trap  lifetime  of  ~30  ms,  and  we  es¬ 
timate  an  atom  temperature  of  ~  30  pK  from  time  of  flight  mea¬ 
surements.  From  the  trap  simulations  (details  are  in  SI  Text),  we 
infer  that  the  atoms  are  confined  to  a  region  145  nm  above  the 
surface  with  dimensions  Axa  ~  ±6  pm  and  A yA  —  Aza  ~  ±30  nm. 
The  simulations  predict  that  more  energetic  atoms  escape  the  trap 
and  collide  into  the  structure,  because  the  weakest  direction  of  the 
trap  is  along  the  diagonals  of  the  y-z  plane  due  to  Casimir- 
Polder  forces. 

To  estimate  the  average  number  of  trapped  atoms,  we  measure 
the  superradiant  atomic  decay  rate  when  the  atom  frequency  vdi  is 
tuned  to  the  first  resonance  v\  of  the  PCW  (Fig.  2C)  (21).  Because 
of  the  strong  dissipative  interactions  between  the  atoms  and  with 
Ad  «  0,  the  collective  decay  rate  is  enhanced  compared  with  the 
single-atom  decay  rate,  and  we  infer  an  average  atom  number  of 
N  =  3.0  ±  0.5  (SI  Text).  In  the  low-density  limit  N  <C  1,  the  mea¬ 
sured  decay  rate  corresponds  to  that  of  a  single  atom.  We  then 
measure  a  GM  decay  rate  Tid  =  (1.5  ±0.2)  To,  which  is  in  good 
agreement  with  the  FDTD  simulations  at  the  calculated  trap  lo¬ 
cation  (SI  Text). 

After  the  atoms  are  loaded  into  the  trap,  we  send  a  weak  5-ms  probe 
beam  Ep  with  frequency  vp  in  either  the  TE  or  TM  GM  through  the 


PCW  and  record  the  transmitted  intensity  \t(vp)  Ep(vp)\2.  The  probe 
beam  scans  near  the  Cs  6Si/2,F  =  3  -»  6Pi/2,F'  =  4  transition. 
Each  experimental  cycle  runs  at  a  fixed  detuning  Aa  =  vp  -  vDi 
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Fig.  3.  Transmission  spectra  of  the  PCW  (A)  without  and  ( B-D )  with 
trapped  atoms.  (/A)  Measured  (black)  and  FDTD-simulated  (blue)  trans¬ 
mission  spectra  of  the  PCW  without  atoms  as  a  function  of  the  probe 
detuning  from  the  band  edge  frequency,  5BE  =vp  -  vBE.  There  is  a  minimum 
extinction  of  25  dB  for  the  transmitted  signal  because  of  fabrication  im¬ 
perfections.  (B-D)  Transmission  spectrum  for  A/  =  3.0±0.5  trapped  atoms 
vs.  probe  detuning  aa  =  i/p-i/di  at  several  frequencies  around  the  band 
edge.  The  solid  lines  are  fits  using  the  transmission  model  in  Eq.  4  aver¬ 
aged  over  atom  positions  and  different  atom  numbers.  In  B,  the  Cs  D ^  line 
is  aligned  to  the  first  cavity  resonance  v1#  resulting  in  symmetric  spectra  for 
both  the  TE  (black;  •)  and  TM  (gray;  A)  modes.  The  TE  spectra  in  C  are  for 
frequencies  on  the  negative  side  (v_;  •)  and  positive  side  (v+;  A)  of  the  v\ 
resonance.  The  TE  spectra  in  D  are  taken  at  the  band  edge  (vBE;  •)  and  60 
GHz  (vBG;  A)  into  the  bandgap.  The  asymmetry  of  the  line  shapes  in  C  and 
D  implies  a  large  ratio  of  coherent  to  dissipative  interactions. 
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relative  to  the  free  space  atomic  transition  frequency  vDi-  We 
observe  little  change  of  signal  during  the  5 -ms  probing  time,  sug¬ 
gesting  that  the  atom  number  is  approximately  constant  over  this 
interval.  The  band  edge  of  the  PCW  is  tuned  thermally  by  shining 
an  external  laser  onto  a  corner  of  the  chip,  where  its  light  is 
absorbed  by  the  Si  substrate.  Hence,  the  Cs  Di  line  can  be  aligned 
to  be  either  outside  or  inside  the  bandgap  with  an  uncertainty 
5v  ~  5  GHz.  The  transmission  for  each  data  point  is  normalized  by 
the  transmission  with  no  atoms  (\toEp  |  ),  resulting  in  a  measure¬ 
ment  of  T/T0  =  \t/t0\2.  The  logarithm  of  the  measured  and  simu¬ 
lated  transmission  spectra  with  no  atoms  T0  =  \to(vp)\  is  shown 
in  Fig.  3 A. 

Examples  of  transmission  spectra  with  atoms  are  shown  in  Fig. 
3  B-D.  Note  that  the  spectra  are  shifted  12.5  MHz  because  of 
both  the  alternating  current  (AC)  Stark  shift  of  the  dipole  trap 
and  the  modified  Lamb  shift  induced  by  the  non-GMs  of  the 
PCW.  Notably,  the  transmission  spectra  at  the  first  cavity  reso¬ 
nance  v\  exhibit  a  characteristic  Lorentzian  “dip,”  and  they  be¬ 
come  more  asymmetric  as  the  frequency  moves  into  the  bandgap. 

Transmission  Model 

We  have  developed  a  model  to  extract  quantitative  values  for 
collective  decay  rates  and  frequency  shifts  from  these  atomic 
transmission  spectra  (39).  Although  the  formalism  of  waveguide 
(40)  and  CQED  (41)  is  well-suited  for  describing  atoms  coupled 
to  uniform  waveguides  and  cavities,  it  is  not  general  enough  to 
capture  the  rich  physics  of  atomic  interactions  in  the  vicinity  of  a 
PCW.  Instead,  we  describe  our  system  by  using  a  spin  model  in 
terms  of  the  classical  electromagnetic  Green’s  function,  in  which 
the  atoms  (or  “pseudospins”  for  ground  and  excited  states)  in¬ 
teract  through  the  emission  and  reabsorption  of  guided  photons 
(42-44). 

The  electromagnetic  Green’s  tensor  G(r,  r/,  co)  is  related  to  the 
electric  field  E(r,  co)  emitted  by  a  dipole  p,  oscillating  at  frequency  co 
at  position  r \  by  E(r,  co)  =  |i0oo2G(r,  r,,  co)  •  p,  (43,  45).  The  dipole 
moment  operator  for  atom  i  is  decomposed  into  p,  =  didlge  +  d*d^, 
where  d;  is  the  dipole  matrix  element  and  6lge  =  \g)(e\  is  the  atomic 
coherence  operator  between  the  ground  and  excited  states.  The 
spin  model  describes  a  system  of  A  atoms  coupled  to  and  driven  by 
a  GM  of  the  PCW.  In  the  low-saturation  and  steady-state  regime, 
expectation  values  for  the  atomic  coherences  (&ge  =  (&ge))  are  de¬ 
scribed  by  a  linear  system  of  equations  (39)  (SI  Text): 

( ~  r\  ^  a 

(4*  +  i  yj  <5lge  +  gij  GJge  =  —Q,  [1] 

where  Aa  =  2tzAa  =  2n(vp  - vdi)  is  the  detuning  between  the 
probe  and  the  atomic  angular  frequencies,  Q  is  the  classical  drive 
(Rabi  frequency)  for  the  ith  atom  due  to  the  GM  input  field,  and 
gij=J ID  +  ir1D/2,  where  Jl(D  =  \i0(e>j/h  cl*  •  Re  G(r„  r ap)  ■  d,  and 
r^D  =  2| i0  or / h  d*  •  Im  G(ry,  i)-,  o)p )  •  d;.  Each  atom  can  also  decay 
into  non-GMs,  including  free  space,  with  a  decay  rate  T'.  The 
appearance  of  the  real  and  imaginary  parts  of  the  Green’s  function 
in  the  coherent  and  dissipative  terms  has  the  classical  analog  that 
the  in-phase  and  out  of  phase  components  of  a  field  with  re¬ 
spect  to  an  oscillating  dipole  store  time-averaged  energy  and 
perform  time-averaged  work,  respectively.  Because  the  first 
term  in  Eq.  1  is  diagonal,  the  atomic  coherences  can  be  un¬ 
derstood  in  terms  of  the  eigenvalues  for  §  =  {1,  mmm,N}  and 
the  eigenfunctions  of  the  matrix  g,  which  has  elements  that  are 
gij ;  the  real  and  imaginary  parts  of  {^}  correspond  to  frequency 
shifts  and  GM  decay  rates,  respectively,  of  the  collective 
atomic  mode 

The  transmission  spectrum  can  be  expressed  in  terms  of  the 
eigenvalues  of  g  as  (39)  (SI  Text) 


Aa  +  ir'/2  \ 
tofa)  f=\\AA+  ir/2  +  At)’ 

where  ^(A^)  is  the  transmission  without  atoms.  In  the  case  of  a 
single  atom  i,  the  only  eigenvalue  is  proportional  to  the  self- 
Green’s  function,  X^  =gu,  which  implies  that  the  transmission  spec¬ 
trum  is  a  direct  measurement  of  the  self-Green’s  function  at  the 
atom’s  position.  For  noninteracting  atoms,  the  off-diagonal  ele¬ 
ments  of  g  are  zero,  and  thus,  the  eigenvalues  are  single-atom 
quantities,  =gu,  because  there  is  no  cooperative  response. 

In  contrast,  for  interacting  atoms,  the  off-diagonal  elements 
are  nonnegligible,  and  there  is  a  cooperative  response.  In  par¬ 
ticular,  for  the  atomic  frequency  inside  the  bandgap  of  a  pho¬ 
tonic  crystal,  the  elements  g^  are  well-approximated  by  (28) 


where  the  cosine  factors  arise  from  the  Bloch  mode,  and  the 
decay  length  1/k*  is  caused  by  the  exponential  decay  of  the  field 
and  results  in  a  finite  range  of  interaction.  For  an  infinite  pho¬ 
tonic  crystal,  TiD  =  0,  because  the  light  is  localized,  and  there  is 
no  dissipation  through  the  GM.  However,  for  a  finite  PCW  of 
length  L,  the  GM  dissipation  riD~e_KxL  is  finite  because  of 
leakage  of  the  mode  out  of  the  edges  of  the  structure. 

In  the  limit  where  the  interaction  range  1/k*  is  much  larger  than 
the  separation  dxy  =  \xt  —Xj\  of  the  atoms,  kx  dxtj  <  kx  Ax  a  1, 
the  GM  input  field  couples  predominantly  to  a  single  collec¬ 
tive  “bright”  mode  of  the  system  with  eigenvalue  Xb  =  Ya=iSh  = 
J2i=i  (Ad  +  i  r“D/2).  Formally,  when  kx  =  0,  the  matrix  g  is  sepa¬ 
rable  \gij  =  UiUj  with  Ui  <xeos(jDC//fl)]  and  therefore,  only  has  one 
nonzero  eigenvalue.  In  this  single  bright  mode  approximation,  the 
transmission  spectrum  is  given  by 


to{AA) 


A^+iT'/2 

(A/i  +  E4d)  +  i  (r'  +  E  Eid)  j 2 


[4] 


We  have  confirmed  numerically  that  this  single  bright  mode 
picture  is  valid  within  the  limits  of  our  uncertainties  for  the  range 
of  frequencies  of  the  measured  spectra  in  Fig.  3.  In  particular,  at 
the  largest  detuning  into  the  bandgap  ABe  =  60  GHz,  we  have 
kx  Axa  —  0.2.  However,  for  atomic  frequencies  farther  away  from 
the  band  edge,  this  approximation  eventually  breaks  down  (e.g., 
at  the  bandgap  center,  kx  Axa  —  1.5). 

The  single  bright  mode  approximation  is  also  valid  in  conven¬ 
tional  CQED.  The  Green’s  function  matrix  is  then  given  by 
gij  =  (J\ d  +  ir id/ 2)  cos (kcXi ) cos (kcXj ) ,  where  kc  is  the  wavevector  of 
the  standing  wave  cavity.  In  this  case,  J\d  cx  Ac/(1  +  A^/y^)  and 
Tid  ocyc/(l  +  A^/y^),  where  Ac  is  the  detuning  from  the  cavity 
resonance  and  yc  is  the  cavity  linewidth.  Importantly,  the  ratio  of  the 
imaginary  dissipative  coupling  rate  to  the  real  coherent  coupling 
rate  falls  off  with  inverse  detuning,  Rcqed  =  T]  d  jJ\  d  =  Yc  /  Ac  for 
large  Ac,  whereas  in  a  PCW  bandgap,  the  fall  off  is  exponential  with 
detuning  from  the  band  edge. 

Analysis  of  Measured  Spectra 

Eq.  4  provides  a  direct  mapping  between  the  observed  trans¬ 
mission  spectra  in  Fig.  3  B-D  and  the  electromagnetic  Green’s 
function  of  the  PCW.  In  particular,  the  line  shape  is  Lorentzian  for 
purely  dissipative  dynamics  (/“D  =  0).  This  line  shape  is  precisely 
what  occurs  at  the  frequency  of  the  first  cavity  mode  v\,  as  shown 
in  Fig.  3 B.  When  the  GM  band  edge  frequency  is  moved  toward 
the  atomic  resonance  vdi,  the  dispersive  interactions  are  switched 
on,  and  the  transmission  line  shape  becomes  asymmetric,  displaying 


4  Of  6  I  www.pnas.org/cgi/doi/i0.i073/pnas.i6^^|^UT|ON  A:  Distribution  approved  for  public  release. 


Hood  et  al. 


Fig.  4.  (A)  Peak  dissipative  interaction  rate  A/T1D  (green)  and  coherent  rate 
/Wid  (blue)  around  the  band  edge.  With  N  determined  from  independent 
decay  rate  measurements,  the  values  for  r-\DlJ-\D  are  found  from  fits  of  the 
transmission  model  in  Eq.  4  to  the  measured  atomic  spectra  and  normalized 
by  the  free  space  decay  rate  r0  =  27c x 4.56  MHz  for  the  Cs  D<\  line.  The  lines 
are  the  predictions  from  a  numerical  model  based  on  ID  transfer  matrices. 
(B)  The  measured  and  calculated  ratios  1Z  =  r-[D/J^D.  The  average  of  the  two 
points  in  the  bandgap  gives  a  ratio  of  the  dissipative  to  coherent  coupling 
rate  IZ  =  0.05  ±0.1 7.  B,  Inset  is  a  comparison  of  IZ  for  the  PCW  calculation 
(solid  line)  and  CQED  model  (dashed  line).  From  the  measured  linewidth  of 
the  first  cavity  resonance,  yc  =  60±8  GHz,  CQED  predicts  that  7£Cqed  =  yc/Ac, 
where  Ac  =  (vp-v1).  Note  that  -V1D  is  plotted  to  more  readily  compare  r-|D 
and  7-id  as  the  band  edge  is  approached. 

a  Fano-like  resonance  (46),  which  can  be  observed  in  Fig.  3  C  and 
D.  The  appearance  of  an  asymmetry  in  the  atomic  spectra  directly 
reveals  a  significant  coherent  coupling  rate  /id,  which  is  evident  for 
frequencies  that  are  in  the  bandgap  region. 

For  all  relevant  frequencies,  the  spectra  for  the  TM  GM  are 
approximately  symmetric,  since  J™  C  T™  «  T'  for  this  GM 
polarization.  An  example  of  a  TM  spectrum  is  shown  as  the  gray 
curve  in  Fig.  3 B.  Because  the  TM  bandgap  is  so  far  detuned,  the 
TM  spectra  are  insensitive  to  ABe  and  serve  as  a  calibration  signal. 
Using  a  waveguide  transmission  model,  we  fit  the  TM  transmission 
spectra  and  extract  a  TM  GM  decay  rate  of  T™  =  (0.045  ±  .01)  T0. 
This  rate  is  ~  30  times  smaller  than  the  TE  GM  decay  rate  F iD  at 
the  first  resonance  v\.  The  ratio  r]§/r™  «  30  is  explained  well  by 
the  expected  slow-light  and  cavity  enhancement  of  the  PCW  de¬ 
scribed  in  ref.  21  and  SI  Text.  From  the  TM  fits,  we  also  measure 
r'  =  2jrx9.1  MHz,  which  because  of  inhomogeneous  broadening, 
is  larger  than  the  value  T'  =  2jt  x  5.0  MHz  predicted  from  FDTD 
numerical  calculations  (SI  Text).  While  tuning  the  band  edge  to 
move  the  atomic  frequency  edi  into  the  bandgap,  TM  spectra  are 
measured  to  confirm  in  situ  that  the  average  atom  number  is  ap¬ 
proximately  constant  over  the  course  of  the  measurements  of 
TE  spectra. 


To  obtain  quantitative  values  for  the  collective  frequency  shifts 
and  decay  rates  by  fitting  the  TE  atomic  spectra  to  the  spin 
model,  we  must  account  for  the  fluctuations  in  atom  number  and 
position  along  the  x  axis.  As  depicted  in  Figs.  L4  and  2 C,  trapped 
atoms  are  approximately  free  to  move  along  the  axis  of  the  de¬ 
vice  (SI  Text).  Their  coupling  rates  are  thus  modulated  by  the  fast 
oscillation  of  the  Bloch  function,  which  near  the  band  edge,  is 
approximately  given  by  Eq.  3,  r“D (x/)  =  Tid  cos2 (xiit/a),  and 
JiD(xi)  =Jm  cos2(x;jc/fl).  Here,  Tid  and  /id  are  the  peak  values. 
Furthermore,  although  we  know  the  average  atom  number 
TV  =  3.0  ±0.5  atoms  from  independent  decay  rate  measurements  (SI 
Text),  the  atom  number  for  each  experiment  follows  an  unknown 
distribution.  To  model  the  experimental  transmission  spectra,  such 
as  in  Fig.  3,  we  average  the  expression  in  Eq.  4  over  the  atom 
positions  {x;}  along  the  Bloch  function  and  assume  a  Poisson 
distribution  (TV)  for  the  atom  number  TV.  We  extract  peak  values 
Tid  and  /id  and  plot  the  resulting  cooperative  rates  NT iD  and 
TV/id  in  Fig.  4 A.  In  particular,  at  the  first  resonance  v\,  the  fitted 
single-atom  GM  decay  rate  is  TiD  =  (1.4  ±  0.2)  To,  which  is  in  good 
agreement  with  the  decay  time  measurements  Tid  =  (1.5  ±  0.2)  To. 
More  generally,  we  find  good  agreement  between  our  measure¬ 
ments  and  our  model  for  the  transmission,  as  shown  in  Fig.  3. 

The  ratio  IZ  =  Tid  //id  is  shown  in  Fig.  4 B.  Because  of  the  ev¬ 
anescent  nature  of  the  field  in  the  bandgap,  IZ  decays  exponen¬ 
tially  with  increasing  detuning  into  the  bandgap,  IZ~e~KxL,  where 
kx  oc  v/Abe  (28).  As  displayed  in  Fig.  4 B,  Inset,  the  ratio  between 
the  GM  decay  rate  F iD  to  the  GM  frequency  shift  /iD  diminishes 
much  faster  than  would  be  the  case  in  traditional  settings,  such  as 
CQED,  for  which  T^cqed  =  yc/Ac,  where  yc  is  the  cavity  linewidth 
and  A c  is  the  detuning  from  the  cavity  resonance.  Indeed,  by 
performing  an  average  of  the  last  two  measured  frequencies  in  the 
bandgap,  we  obtain  IZ  =  0.05  ±  0.17,  whereas  7^cqed  =  0.30  ±  0.04, 
where  we  have  taken  the  cavity  linewidth  to  be  a  value  consis¬ 
tent  with  the  linewidth  of  the  first  cavity  mode  of  the  PCW 
(yc  =  60  ±  8  GHz).  We  can  then  infer  that  the  ratio  of  dispersive 
to  dissipative  rates  for  GM  atom-atom  interactions  (i.e.,  1/IZ) 
is  significantly  larger  than  is  the  case  in  conventional  optical 
physics  (e.g.,  CQED). 

Beyond  the  detailed  modeling  involving  Eq.  4  averaged  over 
fluctuations  in  atom  number  and  position,  we  also  fit  the  spectra  with 
a  generic  transmission  model  with  no  averaging,  as  shown  in  SI  Text. 
We  find  that  the  effective  values  for  the  GM  decay  rate  and  fre¬ 
quency  shift  are  related  to  NT iD  and  NJ\D  in  Fig.  44  by  a  simple 
scale  factor  related  to  the  averaging  of  the  Bloch  function  cos2(ji x/a). 

Despite  favorable  scaling  between  the  collective  frequency  shifts 
and  the  GM  decay  rates,  there  is  still  one  obstacle  to  overcome 
toward  purely  dispersive  atomic  interactions,  namely  atomic  emis¬ 
sion  into  non-GMs  (characterized  by  T').  For  this  PCW  structure, 
the  FDTD-simulated  value  of  this  decay  rate  is  F'  ~  1.1  r0  (21)  for 
the  relevant  frequencies  of  our  experiment.  Fortunately,  it  has  been 
shown  that  suitable  engineering  of  a  wide  variety  of  nanophotonic 
structures  can  lead  to  significant  reductions  in  F'/Fo  (47).  For  ex¬ 
ample,  ref.  1  reviews  possibilities  to  achieve  F'  ~  O.lTo. 

Concluding  Remarks  and  Outlook 

In  conclusion,  we  report  the  initial  observation  of  cooperative  atom 
interactions  in  the  bandgap  of  a  PCW.  By  tuning  the  band  edge 
frequency  of  the  PCW,  we  are  able  to  modify  the  interactions  be¬ 
tween  the  atoms  that  are  trapped  close  to  the  device,  reducing  the 
dissipative  relative  to  coherent  coupling  for  frequencies  inside  the 
bandgap  of  the  PCW.  Equipped  with  a  theoretical  model  based  on 
the  electromagnetic  Green’s  function  of  the  alligator  PCW,  we  infer 
quantitative  values  for  the  collective  frequency  shifts  and  decay  rates 
experienced  by  the  atoms.  Moreover,  we  infer  a  suppression  of  the 
dissipative  interactions  with  respect  to  the  coherent  ones  several 
times  larger  than  is  customarily  obtained  in  atomic  physics.  This 
measurement  provides  the  first  stepping  stone  toward  the  realization 
of  quantum  many-body  physics  in  bandgap  systems. 
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Moreover,  near-term  extensions  of  our  experiment  open  the  door 
to  exploring  new  physical  scenarios  by  using  atoms  coupled  to  PCWs. 
By  trapping  the  atoms  at  the  center  of  the  device  with  GMs  (47),  we 
expect  a  sixfold  increase  to  both  coupling  strengths  /id  and  FiD  rel¬ 
ative  to  r.  Moreover,  by  probing  the  atoms  with  the  Cs  D2  line  tuned 
to  the  upper  band  edge,  where  the  intensity  at  the  position  of  the 
atoms  is  larger,  we  expect  a  further  improvement  by  a  factor  of  two. 
Combining  these  two  effects,  we  expect  a  significant  enhancement  of 
interactions  through  GMs  compared  with  conventional  free  space 
interactions,  namely /id,  Hd  >  10  x  P.  This  improvement  could  enable 
investigations  of  new  paradigms  for  atom-photon  interactions  (28,  29, 
36),  including  the  recently  proposed  multiphoton  dressed  states  (26, 27). 

Note.  After  the  submission  of  this  manuscript,  ref.  48  reported 
measurements  of  transmission  spectra  for  a  superconducting  qubit 
placed  within  the  bandgap  of  a  microwave  photonic  crystal. 
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SI  Text 

In  our  results  in  the  main  text,  we  measure  collective  frequency  shifts 
and  decay  rates  for  atoms  trapped  near  a  PCW.  In  our  previous  work  in 
ref.  21,  we  trapped  multiple  atoms  in  an  optical  dipole  force  trap 
above  the  PCW.  We  operated  with  the  atomic  frequency  outside  the 
bandgap  in  a  regime  with  large  decay  rate  Tid  and  small  coherent 
coupling  rate  Jyd.  By  varying  the  density  and  observing  the  super¬ 
radiant  decay  of  the  atoms  =^sr(N)  +  r]D  +  F,  we  inferred  the 

single-atom  GM  decay  rate  T iD  and  the  average  number  of  atoms  N. 
Importantly,  this  measured  single-atom  decay  rate  F1D  agreed  well 
with  the  FDTD  simulations  at  the  calculated  trap  location.  This 
good  agreement  is  in  part  because  of  the  nanometer-scale  ac¬ 
curacy  in  which  the  alligator  PCWs  are  fabricated,  which  is  re¬ 
quired  for  both  the  band  edge  alignment  and  the  device  quality. 

In  our  paper,  the  band  edge  of  the  PCW  is  tuned  around  the 
resonance  frequency  of  the  atoms,  and  we  observe  the  dominance  of 
the  GM-coherent  coupling  rates  /id  over  the  dissipative  coupling 
rates  TiD,  which  is  associated  with  atomic  radiative  processes 
for  operation  within  the  bandgap.  To  extract  quantitative  values  for 
these  parameters  from  our  measurements  of  transmission  spectra 
for  atoms  trapped  along  a  PCW,  we  have  developed  theoretical 
techniques  based  on  Green’s  functions  for  the  PCW,  which  are  new 
to  atomic  physics.  As  in  ref.  21,  the  average  number  of  atoms  N  is 
measured  by  way  of  transient  decay.  Our  principal  finding  relates 
to  the  turning  off  of  the  GM  decay  rate  T iD,  which  in  the  bandgap, 
is  predicted  to  be  exponentially  suppressed,  while  nonetheless, 
retaining  appreciable  coherent  processes  described  by  J\ d. 

For  the  spectra  in  our  paper,  the  transmission  through  the 
device  decreases  exponentially  in  the  bandgap,  and  more  time  is 
required  to  measure  the  transmission  spectra  compared  with  our 
work  in  ref.  21.  Unfortunately,  Cs  slowly  coats  the  PCW  during 
the  measurement,  both  degrading  the  device  quality  and  shifting 
the  band  edge  out  of  the  thermal  tuning  range.  As  a  result,  each 
device  only  has  a  limited  lifetime  for  making  transmission  mea¬ 
surements.  For  our  experiment,  we  first  repeated  superradiance 
measurements  outside  the  bandgap  at  the  first  resonance  v\  of  the 
PCW  to  determine  the  average  number  of  atoms  N  and  the  single¬ 
atom  GM  decay  rate  Tid  and  show  that  the  atoms  behave  as 
a  collective  emitter.  Then,  with  an  average  number  of  N  ~  3,  we 
measured  transmission  spectra  as  the  atomic  frequency  is  shifted 
into  the  bandgap.  We  simultaneously  measured  the  TM  spectra  to 
verify  that  the  atom  number  is  constant  over  the  course  of  the 
measurements  of  the  TE  spectra. 

1.  Alligator  PCW  Design  and  Fabrication.  The  schematic  of  the  device 
is  shown  in  Fig.  SL4.  Light  is  coupled  into  and  out  of  the  device 
by  mode-matching  the  output  of  an  optical  fiber  to  that  of  a 
terminated  rectangular-shaped  waveguide  on  both  sides  of  the 
device  (33).  The  fibers  are  glued  permanently  in  etched  v 
grooves  at  optimized  coupling  positions.  The  design  and  fabri¬ 
cation  of  the  alligator  PCW  are  detailed  in  ref.  33.  The  PCW  is 
fabricated  on  a  200-pm  Si  chip  coated  with  a  200-nm-thick  SiN 
film.  The  SiN  device  is  suspended  across  a  2-mm-wide  window 
after  the  Si  substrate  beneath  it  is  removed,  as  shown  in  the  image 
in  Fig.  SI B.  The  window  allows  optical  access  for  the  trapping  and 
cooling  of  atoms  around  the  device. 

The  dielectric  TE  mode  band  edge  (eBe)  is  aligned  to  within 
200  GHz  of  the  Cs  Dx  line  (vdi  =335.12  THz)  by  a  low-power 
inductively  coupled  reactive  ion  CF4  etch.  The  directional  etch  thins 
the  SiN  layer  at  a  rate  of  3  nm/min  until  a  transmission  measure¬ 
ment  confirms  alignment  of  the  band  edge.  The  final  geometric 
dimensions  of  the  device  used  in  the  text  are  given  in  Fig.  SIC. 


For  the  experiment,  the  chip  is  placed  at  the  center  of  an  ul- 
trahigh  vacuum  chamber,  and  the  optical  fibers  exit  through  Teflon 
fiber  feedthroughs.  We  measure  the  transmission  through  a  device 
using  a  superluminescent  diode  as  the  source  and  an  optical 
spectrum  analyzer  as  the  detector.  The  measured  transmission  and 
reflection  spectra  are  shown  in  Fig.  S24.  The  transmission  spectra 
near  the  lower  (dielectric)  and  upper  (air)  band  edges  are  com¬ 
pared  with  an  FDTD  simulation  in  Fig.  S2  B  and  C. 

2.  Alligator  Dispersion  Relation  from  Scattering  Images.  Here,  we  de¬ 
scribe  the  analysis  performed  for  the  PCW  dispersion  relations  in  Fig. 
2 E.  We  send  a  single-frequency  laser  beam  through  the  device  and 
image  the  scattered  light  with  a  microscope.  We  integrate  the  image 
over  the  width  of  the  PCW  to  produce  a  single  plot  of  intensity  vs. 
position.  Then,  we  scan  the  laser  frequency  around  the  lower  band 
edge  to  produce  a  2D  plot  of  scattered  intensity  as  a  function  of 
position  x  along  the  device  and  frequency  v  of  the  input  light. 

The  weak  scattered  light  comes  from  small  fabrication  imper¬ 
fections  or  intrinsic  material  defects  and  serves  as  a  probe  of  the  local 
intensity.  Because  each  scatterer  emits  light  at  a  different  rate,  we 
have  to  normalize  the  scattered  light  by  a  reference  intensity  spec¬ 
trum  in  which  the  intensity  of  the  device  is  known.  For  this  reference 
spectrum,  we  average  over  the  intensities  for  frequencies  far  from  the 
band  edge,  where  the  PCW  behaves  like  a  waveguide  and  the  local 
intensity  in  the  device  is  approximately  constant.  The  normalized 
data  are  shown  in  Fig.  S3,  and  a  zoomed-in  version  is  in  Fig.  24. 

In  the  FDTD  simulation  described  above,  we  calculate  the  in¬ 
tensity  along  the  center  of  the  device  for  frequencies  around  the  band 
edge.  Taking  the  maximum  intensity  in  each  unit  cell  and  normal¬ 
izing  by  the  intensity  in  the  waveguide  regime,  we  produce  Fig.  2 B. 

Next,  we  fit  the  intensity  spectrum  at  a  given  frequency  to  a 
model  to  extract  the  wavevector  for  that  frequency.  Near  the 
band  edge,  the  field  in  an  infinite  PCW  is  well-approximated  by 
E(x)  <xcos(x%/a)eldkxX,  where  dkx  =  n/a-kx  in  the  propagating 
band  ( ABe  <  0)  and  5 kx  =  ikx  inside  the  bandgap  (ABe  >  0)-  The 
edges  of  a  finite  photonic  crystal  reflect  with  Rt  because  of  a  large 
group  index  mismatch  between  the  waveguide  section  and  the 
PCW.  The  resonances  of  the  weak  cavity  result  in  the  cavity-like 
intensity  profiles  seen  at  frequencies  1^2,3, 4, 5  in  Fig.  S3.  The  in¬ 
tensity  at  a  point  x  along  a  finite  photonic  crystal  of  length  L  is 
well-approximated  by  a  model  based  on  the  intensity  in  a  cavity 
with  two  mirrors  of  reflectivity  Rt : 

|£(x)|2=/i  \e^x_RieimxLe-mxx^  [S1] 

where  I\  is  related  to  the  overall  intensity.  This  expression  ignores 
the  fast  oscillations  of  the  Bloch  function,  which  go  as  cos2(xjc /a). 
Note  that  in  the  bandgap  (when  K XL  >  1),  the  intensity  model 
reduces  to  an  exponential  decay:  \E(x)\2  e~2KxX.  Interestingly, 
at  the  band  edge  (5 kx  ->  0  and  Rt->  1),  the  intensity  displays  a 
quadratic  dependence  on  the  position:  \E(x)\2  oc  (L  -x)2. 

For  each  frequency,  we  fit  the  intensity  along  the  nominal  cells 
with  Eq.  SI  and  extract  5 kx.  This  procedure  allows  us  to  map  out 
the  dispersion  relation  5/c*(ABe),  which  we  show  in  Fig.  2 E  for  the 
measured  and  simulated  data.  From  the  simulated  fits,  we  find 
that  the  effective  length  of  the  cavity  is  162  cells,  which  is  slightly 
longer  than  the  150  nominal  cells  and  is  expected  due  to  the 
leakage  of  the  cavity  field  into  the  tapering  sections.  We  use  this 
length  for  the  fits  of  the  measured  data.  Examples  of  the  mea¬ 
sured  and  simulated  intensities  are  shown  in  Fig.  S4.  The  fluctu¬ 
ation  of  the  intensity,  even  after  the  normalization,  is  most  likely 
caused  by  the  spatial  profile  of  Bloch  mode.  The  normalization 
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trace  is  taken  by  averaging  data  for  excitation  frequencies  farther 
away  from  the  band  edge  where  the  Bloch  mode  contrast  is  re¬ 
duced,  whereas  the  data  closer  to  the  band  edge  have  a  large 
Bloch  mode  fringe  visibility.  However,  the  fluctuations  do  not 
affect  the  statistical  fits  at  the  level  of  accuracy  required  for  the 
dispersion  relation  in  this  work. 

The  frequency  for  which  5 kx  =  0  is  defined  as  the  band  edge 
frequency  ebe-  To  extract  this  frequency  and  the  curvature  of  the 
dispersion  relation  near  the  band  edge,  we  fit  the  measured  and 
simulated  dispersion  relations  with  a  dispersion  model  (21), 


bkx{y) 


2jc  /(e BE2  -  v)  (EBE  -  v) 

a  V  "  (vbe2  “  vbe)2 


[S2] 


where  ebe  (ebe2)  is  the  lower  (upper)  band  edge  frequency,  and  C, 
is  a  frequency  related  to  the  curvature  of  the  band  near  the  band 
edge.  From  the  measured  data  fits,  the  distance  between  the 
first  resonance  and  band  edge  is  ebe  -  v\  =  133  ±  9  GHz  and 
£  =  227  ±  3  THz.  The  simulated  data  give  ebe  —  v\  =  135.0  GHz, 
and  the  curvature  parameter  is  £  =  226.0  THz.  These  values  are 
in  good  agreement  with  the  dispersion  relation  from  the  eigen¬ 
mode  simulation  of  the  infinite  PCW  in  Fig.  1C,  which  gives 
l,  =  229.1  THz. 

3.  SI  Trap.  In  Fig.  S5/4,  we  show  a  schematic  of  the  SI  trap.  The  SI 
beam  is  nearly  perpendicular  to  the  axis  of  the  device,  has  a 
50-pm  diameter,  and  has  a  polarization  aligned  to  the  axis  of  the 
device  (Fig.  S5A).  The  orange  areas  in  Fig.  S5/4  represent  the 
approximate  localization  of  the  atoms  along  x,y.  By  time  of  flight 
measurements  of  atoms  in  the  dipole  trap,  we  estimate  an 
atomic  temperature  of  ~30  pK.  From  the  beam  waist  and  atom 
temperature,  we  can  infer  that  the  atoms  are  localized  to 
2 Axa  =  12  pm  along  the  x  axis. 

Simulations  of  the  SI  trap  potential  are  shown  in  Fig.  S5  B-D. 
The  simulations  are  performed  for  the  infinite  structure  with 
COMSOL.  The  trap  depth  is  calibrated  with  the  12-MHz  AC 
Stark  shift  measured  from  the  atomic  spectra.  Fig.  S 5B  shows  the 
trap  potential  in  the  y-z  plane.  Atoms  that  are  significantly 
hotter  than  ~  100  pK  are  expected  to  crash  into  the  device  along 
the  diagonal  directions  because  of  Casimir-Polder  forces.  Fig. 
S5C  shows  the  trapping  potential  along  the  z  axis.  Atoms  are 
trapped  at  z  =  240  nm.  Fig.  S 50  shows  the  trap  along  the  x  axis. 
Because  of  the  photonic  crystal,  the  trap  modulates  by  ~  10  pK 
along  the  x  axis,  which  is  significantly  smaller  than  the  estimated 
trap  temperature. 

In  addition  to  the  results  in  Fig.  S5,  we  have  also  carried  out 
numerical  modeling  of  the  optical  trap  using  Lumerical  simula¬ 
tions  (37)  of  the  actual  finite  length  PCW  and  tapers  shown  in 
Fig.  SI.  We  have  as  well  included  Casimir-Polder  potentials  as  in 
ref.  47.  More  details  of  the  trap  are  discussed  in  ref.  21. 


4.  Transmission  Model  and  Atomic  Spectra  Fits.  Here,  we  give  a  more 
detailed  description  of  the  transmission  model  in  the  text,  which 
follows  the  derivation  given  in  ref.  39.  A  system  of  A  atoms  coupled  to 
a  radiation  field  can  be  described  using  formalism  based  on  the 
classical  Green’s  function  (42,  43).  In  the  Markovian  limit,  the  field 
can  be  eliminated  to  obtain  a  master  equation  that  describes  the 
interactions  between  the  atoms:  pA  =  - i/h[H ,  pA\  +  C[pA].  Here,  the 
Hamiltonian  H  gives  the  coherent  evolution  of  the  system: 


N  N 

H  —  h  ^  ^  A/±ole  fl  ^  ^  J |  \^^lg^ge  k 
7=1  74=1 


i= 1 

[S3] 


and  the  Lindblad  operator  C[pA]  gives  the  dissipation  of  the  system: 


^]=E^ 

74=1 

X  (2 d}ge  pA  deg  —  dig  dge  pA  —  pA  dig  dge  ) . 


[S4] 


The  Hamiltonian  and  the  Lindblad  are  expressed  in  terms  of  the 
atomic  coherence  operator  &  =  \g)(e\  between  the  ground  and 
excited  states  of  atom  ;.  The  Hamiltonian  contains  terms  for  the 
free-atom  evolution,  the  coherent  atom-atom  interactions,  and 
the  classical  drive,  respectively;  Aa  =  2kAa  =  2n(vp  -  edi)  is  the 
detuning  between  the  probe  and  the  atomic  angular  frequencies, 
and  £lj  is  the  Rabi  frequency  for  atom  ;  caused  by  the  GM  field. 
The  atom-atom  spin  exchange  rate  /^D  is  expressed  in  terms  of 
the  real  part  of  the  GM  Green’s  function  as 


jji  _ 
JYD  ~ 


•Re  G(r,,r„Mp)'d,, 


[S5] 


where  ®p  =  2nvp,  and  d7  is  the  dipole  matrix  element  of  atom  ;. 
The  Lindblad  term  is  responsible  for  the  dissipative  interactions 
in  the  system,  which  include  atomic  decay  into  non-GMs  (T')  and 
GMs  (I^D).  The  decay  rate  into  the  GM  is  written  in  terms  of  the 
imaginary  part  of  the  Green’s  function  as 


r to  =  2  ^44  j  d* .  Im  G(ry,  r„  o>„)  •  d„  [S6] 

For  low  atomic  density  along  the  PCW,  the  nonguided  emission 
rate  T'  is  not  cooperative  and  is  described  here  as  a  single-atom 
effect,  with  dp  as  the  Kronecker  delta. 

In  the  low  saturation  regime,  the  Heisenberg  equations  for  the 
expectation  value  of  the  atomic  coherences  ((deg)  =ceg)  can  be 
solved  for  with  the  master  equation  leading  to 


n 


where  the  complex  coupling  rate  is 

gij  =4d  +  ^  =  d*  •  G(r„  Tj,  (Dp) •  dj,  [SSI 

which  is  the  Green’s  function  between  atoms  i  and  j  projected  onto 
the  respective  dipole  matrix  elements.  In  the  steady-state  solution, 
the  time  derivative  is  set  to  zero,  and  the  result  is  the  linear  system 
of  equations  for  the  atomic  coherences  given  in  the  text. 

The  electric  field  in  the  system  can  be  expressed  in  terms  of  the 
input  probe  field  E+(r,  odp)  and  solutions  for  the  atomic  coherences 
(39): 


+  i  +  i  ^  ^jSji  [S7] 


<*ge  : 


M  i 


N 

E+  (r,  (Dp)  =  E+  (r,  (Dp)  +  p0coj  ^  G(r,  ry,  (Dp)  •  .  [S9] 

7=1 

An  expression  for  the  transmission  through  a  quasi- ID  structure 
can  be  derived  by  solving  the  steady-state  system  of  equations  in 
Eq.  S7  for  the  atomic  coherences  dge  and  substituting  them  into 
Eq.  S9.  The  expression  can  then  be  simplified  in  the  case  where 
the  dipole  moments  are  real,  in  which  case  $  is  a  complex  sym¬ 
metric  matrix  with  eigenvectors  and  eigenvalues  u^,  and 

where  the  Green’s  function  is  well-represented  by  a  ID  Green’s 
function.  The  final  result  is  (37) 
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t(LA,N)_£rf  A^+ir'/2  \ 
to(  Aa)  \}\^A+i^/2  +  X%), 


[S10] 


where  *o( Aa)  is  the  transmission  without  atoms. 

In  the  bandgap,  the  matrix  g  of  elements  is  well-approximated  by 


SH  = 


[Sll] 


As  discussed  in  the  text,  when  the  interaction  range  1  /k*  is  much 
larger  than  the  separation  distance  (kx| xi  —Xj\  <C  1),  there  is  only 
a  single  atomic  bright  mode  for  which  the  frequency  shift  and 
GM  decay  rate  are  given  by  YaLi^id  and  The  trans¬ 

mission  spectrum  for  N  atoms  in  the  single  bright  mode  approx¬ 
imation  is  given  by 


T(Aa,N)  =  T0(Aa) 


AA  +  iT'/2 


Aa  + 


ir'/2+E(/fc+  fllo/2) 


[S12] 


where  AA=2nAA=2n(vp -vm)  is  the  detuning  between  the 
pump  and  the  atomic  frequency,  and  Tq(Aa)  is  the  device  trans¬ 
mission  when  no  atoms  are  present. 

Explicitly  accounting  for  the  atoms’  positions  by  substituting 
Eq.  Sll  into  Eq.  S12,  we  find  that  the  transmission  is  given  by 


T(Aa,N;xu  . . .  ,xn)/T0(Aa) 


Assuming  N  =  3.0,  which  is  obtained  from  the  atom  decay  rate 
measurement,  we  fit  the  TE  atomic  spectra  with  Eq.  S14  and 
extract  Tid,  A d,  F,  and  Ao  for  each  frequency.  We  show  the 
values  of  Tid  and  Ad  in  Fig.  44.  We  show  the  AC  Stark  shift  and 
nonguided  decay  rate  in  Fig.  S6. 

The  average  of  the  nonguided  decay  rate  T'  for  the  TE  data 
outside  the  bandgap  is  F  =  2jcx9.1  MHz,  which  is  significantly 
larger  than  the  expected  value  from  the  FDTD  simulation, 
F  =  2jcx5.0  MHz.  This  additional  inhomogeneous  broadening 
could  be  caused  by  finite  temperature  of  the  trapped  atoms, 
vector  shifts  from  circular  light  in  the  SI  beam,  atom  density- 
dependent  collisional  broadening,  stray  magnetic  fields,  and 
electric  fields  from  charges  in  the  dielectric.  We  estimate  the 
contributions  individually  and  find  that  they  likely  do  not  explain 
the  extraneous  broadening.  We  note  that  the  estimate  of  tem¬ 
perature  of  trapped  atoms  could  be  improved  in  the  future  (49), 
and  it  may  help  shed  light  on  our  excess  broadening. 

Interestingly,  the  fitted  F  increases  in  the  bandgap  and  is  as 
high  as  F  =  2n  x  16  MHz  for  the  last  measured  point.  One  pos¬ 
sible  explanation  is  that  this  is  because  of  the  breakdown  of  the 
single  bright  mode  approximation,  because  coupling  to  multiple 
collective  atomic  modes  should  result  in  a  broadened  linewidth. 
Another  possibility  is  that,  because  there  is  a  large  extinction  of 
the  TE  mode  in  the  bandgap,  there  might  be  some  mixing  be¬ 
tween  the  TE  and  TM  modes. 

We  also  measure  transmission  spectra  for  the  TM  mode,  which 
have  band  edges  that  are  far-detuned  from  the  Cs  transitions. 
The  transmission  in  this  waveguide  regime  is  described  by  an 
optical  depth  (OD)  model 


A4+ir/2 

2 

[S13] 

II 

8 

-OD 

Aa  +iF/2+  (jiD  +  iTiD/2^cos2 

[‘♦(r-ir)  J 

[S15] 


We  have  defined  KA  =  Aa  +  Ao  to  account  for  the  AC  Stark  shift 
Ao  of  the  atoms  because  of  the  dipole  trap. 

To  accurately  model  the  experimental  conditions,  we  average 
the  transmission  model  over  atom  positions  and  atom  number. 
During  a  single  measurement,  the  atoms  are  free  to  move  along 
the  length  of  the  device  over  the  range  2Axa  as  in  Fig.  S5A,  evenly 
sampling  the  Bloch  function.  We  let  (T(Aa,N;xi,  . .  -,xN))x  be 
an  average  over  all  positions,  that  is, 


(T(Aa,N;xu  ...,xn)\ 

a 

=  T0(AA)Jdx1...dxJ 


A 

a  +  iF'/2 

A>iT72+El 

j= i 

(Ad  +  iTio/2^ 

|  cos2 1 

(?) 

We  repeat  the  measurement  multiple  times  for  each  frequency 
Aa.  Each  experiment  can  have  a  different  number  of  atoms,  and 
therefore,  we  average  the  transmission  expression  over  a  Poisson 
distribution  which  is  a  function  of  the  average  atom 

number  N.  The  transmission  model  averaged  over  both  atom 
positions  and  atom  numbers  is  given  by 

(T(Aa,N;xu  ...,xN))xp 

=  T(s(Aa)  Y,  Pm(N)(T(Aa,N-,Xi,  . . .  ,xN))x.  [S14] 

N 

This  expression  is  the  final  form  of  the  transmission  model  that  we 
use  to  fit  the  atomic  spectra. 


where  the  resonant  OD  is  given  by  OD  =  2ATT™/F.  We  fit  the 
TM  spectra  with  this  model  and  extract  F,  Ao,  and  T™  (as¬ 
suming  N  =  3).  The  values  of  F  and  Ao  are  shown  with  the  TE 
data  in  Fig.  S6.  The  averaged  T™  value  is  0.044  To,  which  is 
~30  times  smaller  than  Fd  for  the  TE  mode  at  the  first  res¬ 
onance  v\  and  clearly  shows  the  enhanced  interaction  because 
of  the  PCW. 


5.  Simple  Transmission  Model.  In  the  text,  we  fit  atomic  transmission 
spectra  with  the  averaged  transmission  model  from  Eq.  S14  to 
extract  the  peak  GM  decay  rate  Tid  and  frequency  shift  Ad-  In 
this  section,  we  fit  the  spectra  with  a  transmission  model  that 
involves  no  averaging,  and  we  extract  an  effective  decay  rate  T®^ 
and  frequency  shift  which  will  be  smaller  than  the  corre¬ 
sponding  peak  values  because  of  the  averaging  of  the  cos2(j xx/a) 
Bloch  function  as  the  atoms  move  along  the  x  axis  of  the  trap.  In 
the  single  bright  mode  approximation  discussed  in  the  text,  the 
transmission  for  a  single  collective  mode  with  total  decay  rate  A 
and  frequency  shift  B  is  given  by 


T(Aa)  =  A^+iF/2 

Tq(Aa)  Aa  +B  +  i  (F  +A)/2 


[S16] 


Here,  the  detuning  KA  includes  the  AC  stark  shift  KA  =  Aa  +  Ao. 
Because  the  average  number  of  atoms  N  «  3  is  measured  inde¬ 
pendently  in  a  decay  rate  measurement,  the  collective  rates  A 
and  B  are  related  to  the  effective  rates  by  A  =AT^  and 
B=NJf^.  Examples  of  the  fitted  spectra  for  atoms  outside  and 
inside  the  band  gap  are  shown  in  Fig.  S7.  The  translucent  lines  in 
Fig.  S7  are  the  expected  signals  for  average  atom  numbers  of 
N  =  1  and  N  =  9. 
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The  fitted  values  of  A  and  B  are  plotted  for  each  detuning 
from  the  band  edge  Abe  in  Fig.  S&4.  The  results  are  qualitatively 
similar  to  the  corresponding  plot  in  Fig.  44,  except  that  the  ef¬ 
fective  rates  A  =A/T^  and  B  =TVT®d  are  scaled  down  by  r|  =  0.42 
because  of  the  modulation  of  the  Bloch  function  cos2(j vc/a).  The 
solid  lines  in  Fig.  S&4  are  the  same  theoretical  curves  as  in  Fig. 
44,  except  that  they  are  scaled  by  r|  =  0.42. 

The  ratio  oiA/B  =  T^/J^  is  plotted  in  Fig.  S 8B.  Because  the 
scale  factors  rj  cancel,  the  result  is  in  good  agreement  with  the 
corresponding  plot  of  7£  =  FiD//id  in  Fig.  4 B.  The  black  theory 
curve  in  Fig.  S8 B  is  the  same  as  in  Fig.  4 B.  Whereas  the  peak 
decay  rate  and  frequency  shift  are  sensitive  to  the  specific  model, 
the  ratio  of  dissipative  to  coherent  coupling  is  mostly  model 
insensitive. 


rates,  we  deduce  that  the  apparent  single-atom  decay  rate  is 
fio  =  (1.12  ±0.14)r'. 

Because  the  atoms  are  randomly  distributed  along  the  x  direction 
in  the  trap,  the  observed  decay  curves  are  results  after  spatially  av¬ 
eraging  the  coupling  rates  riD(x).  Assuming  a  uniform  distribution 
of  TV  atoms  around  the  center  of  the  PCW,  a  more  detailed  model 
specifies  the  form  of  fluorescence  intensity  decay  as  (21) 


6.  Atom  Decay  Measurement.  We  exploit  the  superradiance  of  atoms 
trapped  near  the  alligator  PCW  to  determine  the  mean  atom  number 
TV  and  the  peak  atom  decay  rate  TiD  (at  v\)  into  the  GMs. 

As  established  in  ref.  21,  the  total  exponential  decay  rate  of  the 
atoms  is  rtot  (TV)  =  FSr  (TV)  +  Ttot ,  where  rSR  is  the  TV-dependent 
superradiance  decay  rate,  and  T^j  is  the  observed  single-atom 
decay  rate.  We  note  that,  when TV  <C  1,  rtot  ~  f\J  =  F\D  +  T',  because 
only  the  single-atom  decay  rates  T\D  into  the  GM  and  F'  into  the 
environment  remain;  F'  is  numerically  calculated  to  be  2jcx5.0 
MHz  for  the  Cs  Dx  line  at  the  trapping  site  near  the  PCW  (21). 

We  excite  the  atoms  with  a  weak  resonant  light  pulse  through 
the  GM,  whereas  the  first  resonance  v\  near  the  band  edge  is 
aligned  with  the  Cs  D  j  line.  Pulse  properties  are  as  in  ref.  21.  The 
subsequent  fluorescence  decay  rates  rtot  are  determined  through 
exponential  fits.  By  varying  the  trap  holding  time  tm  after  loading, 
the  mean  atom  numbers  for  the  decay  measurements  are  varied. 
The  decay  rates  are  empirically  fitted  in  an  exponential  fomi  as 
a  function  of  holding  time  tm  (21)  :  r,o,(fm)  =  rSRe  'm/TSR  +r$,  as 
shown  in  Fig.  S9.  From  the  fitted  asymptotic  value  of  the  decay 


where  y  =  Tid/2,  and  4  is  the  modified  Bessel  function.  Numerically 
simulating  the  decay  of  single  atoms  in  the  trap  by  using  1\  (f),  we 
compare  between  the  exponentially  fitted  value  F id  and  the  value  of 
T id  used  for  1\  (t),  which  yields  a  ratio  of  F iD /F id  =  0.81.  This  ratio 
is  consistent  with  the  ratio  of  0.8  ±  0.3  from  measurement  at  long 
hold  time  tm  =  94  ms,  when  single-atom  decay  predominates  (shown 
as  the  asymptote  in  Fig.  S9).  Based  on  the  values  of  F id  deduced 
above,  we  conclude  that  TiD  =  (1.4  ±  0.2)T'. 

At  early  holding  times,  the  atom  number  TV  noticeably  fluctuates 
around  some  mean  values  TV  >  1.  To  capture  this  TV-dependent 
variation,  we  fit  the  decay  curves  by  averaging  In  it)  with  weight 
function  of  Poisson  distribution  probability  (TV)  (21).  The  fitting 
parameter  here  is  TV,  whereas  we  fix  the  value  of  TiD  in  Eq.  S17. 
The  fit  is  consistent  with  TV  =  3.0  ±  0.5  at  tm  =  4  ms  when  we  carry 
out  the  transmission  spectra  measurement.  Based  on  the  trap  life¬ 
time  t  =  30  ms,  we  further  deduce  that  TV  ~  0.1  at  tm  =  94  ms. 

The  linear  TV  dependence  of  superradiance  is  given  by 
Tsr  =  q  TV  •  Tid,  where  r|  =  0.36  ±  0.06  is  some  linear  coefficient 
that  has  a  value  is  consistent  with  that  reported  in  ref.  21. 


Fig.  SI.  Alligator  PCW  chip  and  device  overview.  (A)  Schematic  of  the  entire  device.  The  alligator  PCW  is  at  the  center.  Optical  fibers  (green)  on  both  ends 
couple  light  into  and  out  of  the  waveguide.  The  waveguide  is  surrounded  by  supporting  and  cooling  structures.  ( B )  Image  of  a  10  x  10-mm  PCW  chip.  Multiple 
waveguides  stretch  across  the  window  of  the  chip,  with  the  PCWs  at  the  center  of  the  window.  The  window  provides  optical  access  for  trapping  and  cooling 
atoms  around  the  device.  Reproduced  from  ref.  33,  with  the  permission  of  AIP  Publishing.  (C)  Overview  of  device  variables.  The  lattice  constant  for  the  entire 
device  is  a  =  370  nm.  The  device  dimensions  are  measured  with  an  SEM  and  calibrated  to  the  lattice  constant.  The  device  dimensions  are  w  =  310±10  nm, 
2A  =  262 ±  10  nm,  g  =  220±  10  nm,  Wjnitia|  =  268±  15  nm,  and  ginitja|  =  1 65 ±  10  nm.  The  thickness  of  the  SiN  is  1 85 ±  5  nm.  The  index  of  refraction  for  Si3  N4  is 
n  =  2.0  around  our  frequencies  of  interest. 
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Fig.  S2.  Measured  and  simulated  transmission  and  reflection  spectra.  (/A)  Transmission  (black)  and  reflection  (blue)  spectra  through  the  entire  chip  for  the  TE 
mode  (polarization  in  the  plane  of  the  device).  The  red  dashed  lines  are  the  Cs  D -i  (335.1  THz)  and  D2  (351 .7  THz)  lines.  The  TE  transmission  efficiency  through 
the  entire  device  near  the  dielectric  band  edge  is  ~  23%,  indicating  that  the  single-pass  efficiency  from  the  fiber  to  device  is  approximately  49%.  Most  of  the 
loss  is  caused  by  the  waveguide  to  fiber  coupling  section.  The  gray  line  is  the  TM  transmission  (polarization  perpendicular  to  the  plane  of  the  device).  Note  that 
the  lower  band  edge  of  the  TM  mode  is  visible  at  around  365  THz  but  far  detuned  from  both  Cs  D12  transitions.  ( B  and  C)  TE  transmission  data  are  normalized 
and  compared  with  an  FDTD  simulation  (37).  The  simulation  uses  the  measured  device  parameters  in  Fig.  SI  that  are  adjusted  within  the  uncertainty  of  the 
measurements  so  that  the  positions  of  the  first  resonances  match  those  in  the  measured  spectra. 


Fig.  S3.  Normalized  magnitude  of  the  scattered  electric  field  of  the  PCW  for  frequencies  ABe  =  vp  -vBe  around  the  band  edge.  The  schematic  in  Left  shows  the 
PCW  with  the  number  of  unit  cells  reduced  by  five. 
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Fig.  S4.  The  electric  field  magnitude  in  (A)  the  PCW  at  the  first  resonance  v-\  and  ( B )  the  bandgap  vBG  =  vBE  +  60  GHz.  The  points  show  measured  data,  and  the 
black  lines  are  from  an  FDTD  simulation.  The  electric  field  magnitude  |£|  is  normalized  by  the  electric  field  magnitude  far  from  the  band  edge;  thus,  these  plots 
give  the  enhancement  of  |£|  relative  to  the  waveguide  regime. 


Fig.  S5.  (>4)  Schematic  of  the  atoms  in  the  SI  trap.  Given  the  estimated  atom  temperature  of  30  pK,  we  infer  that  the  atoms  are  confined  to  a  length  of 
2A xA  =  12  pm  along  the  x  axis.  ( B-D )  Far-off-resonance  optical  trap  (FORT)  potentials  for  the  SI  trap  simulation  ( B )  in  the  y-z  plane  (21),  (C)  along  the  z  axis,  and 
(D)  along  the  x  axis. 
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Fig.  S6.  Fitted  values  from  the  averaged  transmission  model  for  TE  (black  circles)  and  TM  (gray  triangles)  spectra.  04)  Fitted  AC  Stark  shift  A0.  (£)  Fitted  P. 


Hood  et  al.  www.pnas.org/cgi/content^hort/i603^^|BUT|ON  A:  Distribution  approved  for  public  release. 


6  of  7 


Fig.  S7.  Fits  of  transmission  spectra  with  the  model  in  Eq.  S16  for  when  the  atomic  resonance  frequency  is  aligned  to  (/\)  the  first  resonance  and  (fi)  in  the 
bandgap.  From  the  decay  rate  measurement,  the  average  number  of  atoms  is  N&3  for  the  full  (central)  curves  in  A  and  B,  while  the  translucent  curves  give  the 
expected  spectra  for  N  =  1  and  N  =  9  atoms. 


Fig.  S8.  (/A)  Fitted  values  for  the  effective  collective  decay  rates  A  and  frequency  shifts  B  for  various  detunings  from  the  band  edge  ABe-  The  solid  lines  are  the 
expected  results  for  the  peak  values  as  in  Fig.  4A  except  scaled  down  by  r|  =  0.42.  (B)  Ratio  A/B  =  T^/J^  along  with  the  theoretical  prediction  for  the  peak 
ratio  r1D/71D  from  Fig.  4 B. 


Fig.  S9.  Total  decay  rates  as  a  function  of  holding  time  tm.  The  red  solid  curve  is  the  empirical  fit,  and  the  dash-dot  line  represents  the  fitted  asymptotic  total 
decay  rate  at  very  long  times.  The  blue  dashed  lines  specify  fitted  error  boundaries.  The  fit  yields  xSR  =  16  ms,  rSR  =  1.5F,  and  the  asymptote  F^/r'  =  2.12±0.14. 
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